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Internal transport barriers (ITB's) observed in tokamaks are described by a purely 
magnetic approach. Magnetic line motion in toroidal geometry with broken magnetic 
surfaces is studied from a previously derived Hamiltonian map in situation of incomplete 
chaos. This appears to reproduce in a realistic way the main features of a tokamak, 
for a given safety factor profile and in terms of a single parameter L representing 
the amplitude of the magnetic perturbation. New results are given concerning the 
Shafranov shift as function of L. The phase space (ip,8) of the "tokamap" describes 
the poloidal section of the line trajectories, where ip is the toroidal flux labelling the 
surfaces. For small values of L, closed magnetic surfaces exist (KAM tori) and island 
chains begin to appear on rational surfaces for higher values of L, with chaotic zones 
around hyperbolic points, as expected. Island remnants persist in the chaotic domain 
for all relevant values of L at the main rational q- values. 

Single trajectories of magnetic line motion indicate the persistence of a central pro- 
tected plasma core, surrounded by a chaotic shell enclosed in a double-sided transport 
barrier : the latter is identified as being composed of two Cantori located on two suc- 
cessive "most-noble" numbers values of the perturbed safety factor, and forming an 
internal transport barrier (ITB). Magnetic lines which succeed to escape across this 
barrier begin to wander in a wide chaotic sea extending up to a very robust barrier 
(as long as L ^ 1) which is identified mathematically as a robust KAM surface at the 
plasma edge. In this case the motion is shown to be intermittent, with long stages of 
pseudo-trapping in the chaotic shell, or of sticking around island remnants, as expected 
for a continuous time random walk. 

For values of L y 1, above the escape threshold, most magnetic lines succeed to 
escape out of the external barrier which has become a permeable Cantorus. Statistical 
analysis of a large number of trajectories, representing the evolution of a bunch of 
magnetic lines, indicate that the flux variable if) asymptotically grows in a diffusive 
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manner as (L 2 t) with a L 2 scaling as expected, but that the average radial position 
r m (t) asymptotically grows as (L 2 t) 1/4 while the mean square displacement around this 
average radius asymptotically grows in a subdiffusive manner as (L 2 t) l/2 . This result 
shows the slower dispersion in the present incomplete chaotic regime, which is different 
from the usual quasilinear diffusion in completely chaotic situations. For physical times 
of the order of the escape time defined by Xm(t v ) ~ 1, the motion appears to be 
superdiffusive, however, but less dangerous than the generally admitted quasi-linear 
diffusion. The orders of magnitude of the relevant times in Tore Supra are finally 
discussed. 

PACS numbers: 52.55.Fa, 05.45. +b, 52.25.Gj, 05.40. +j, 52.35.Fp 
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1 Analytic solution for the Shafranov shift: position ip of the magnetic axis as 
function of K = 2irL in the case w=_L For L = 4.875/27T this formula yields 

i[j(L — 4.875/27r) = 0.07492 . Here r — \f4> represents the radial coordinate, 
with r — \J2 on the edge. 

2 Tokamap phase portrait for L = 0.5/27T with q(0) < 1 (w = 1.05 and N — 
8190). The polar axis (0,0) is the main magnetic axis of the plasma core. 
The position of the elliptic point in the m = 1 island is ip = 0.040246, 8 = 0.5 
and the hyperbolic point is located in ip = 0.024481, 6 = 

3 Regular KAM magnetic surfaces with a Shafranov shift for a small value oi 
the stochasticity parameter L = 0.08, with w = 1. The angle 9 is represented 
modulo 2n 

4 Regular KAM tori for L = 0.4, w = 1, and main rational island chains] 
corresponding to q — 5/4, 4/3, 3/2, 2, 3 and 4 



5 With the same initial conditions as in Fig. 4 (w = 1) with a higher value of the 



stochasticity parameter L = 4.875/27T ~ 0.776, regular KAM tori are observed 



in the central part, while the main rational island chains corresponding to q 



3/2, 2, 3 and 4 are stochastized 23 



3 Partial vue of the phase portrait with six trajectories for L = 4.875/2-7T 



0.776, which shows the structures of rational chains of island remnants and 



the stochastic web expanding from the hyperbolic points 24 

7 Part of the skeleton of island remnants in the phase portait (r = \/2ip _ x \J~2 



as function of 9) for L = 4.875/27T ~ 0-776 : rational chains can be seen 



corresponding to q = 4, 3, 2, 3/2, 4/3, 8/7, 9/8, 12/11 J 25 



B Several trajectories for L = 5.5/27T ~ 0.875.. , yielding a central protected 



plasma core, surrounded by chains of regular islands and KAM surfaces, and 



outside, a stochastized zone allowing the line to escape from this central 



chaotic shell 26 



Phase portrait of the plasma edge, after complete destruction of the KAM 



surfaces. Several trajectories are chosen to represent the details of the four 



island remnants q = m — 4 at a value L = 6/2-7T ~ 0.954... beyond the escape 



threshold. The plasma edge has become permeable and is strongly deformed 



as compared to the unperturbed circle ip = 1. One remarks the presence oi 



four satellites or daughter islands around each main island 27 



10 Single trajectory at L = 5.5/27T ~ 0.875 followed during N = 16380 iterations 



After a long trapping period in the chaotic sea, this trajectory finally escape: 



out of the plasma edge. The initial conditions are r = \/2ip — 0.4 and 9 — 0.025.. 28 



11 Phase portrait for L — 4.875/27T ~ 0.776 as drawn in coordinates 9 and 
= y/2ip by a single trajectory followed during 2.10 9 iterations, with one 
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point represented every 10009 iterations. The inner shell appears slightly 



darker, surrounding the plasma core (white) and the magnetic axis. Long 



ticking occurs here around the five small q = 5/2 islands (dark lines in the 



upper part of the chaotic sea) 29 



12 Detail of a small structure of the phase portrait around one of the five q = 5/2 



islands, where a long sticking stage occurs during 10 8 iterations along the 



trajectory represented in Figs. 11 and 12 30 

13 Radial position of a single trajectory followed during 2.10 9 iterations, as func- 



tion of time for L = 4.875/2?r ~ 0.776,with initial condition ijj = 0.0136125 



9 = 0.033 : different stages of trapping in different zones indicate an inter- 



mittent behaviour. Note that trapping stages are very long at may last for at 



least several 10 7 iterations for this value of the stochasticity parameter, ncai 



the escape threshold. (Such long times near the critical threshold for large 



scale chaotic motion could be reminiscent of critical phenomena 32 



14 Structure of rational surfaces (islands) and chaotic zones below the edge 



KAM torus iV(4, 2) (boundary circle) and the three convergents from bc- 



low to iV(4,2): q = 13/3, q = 35/8, q = 92/21. This detailed drawing foi 



L = 4.875/27T ~ 0.776 represents short iteration times (N = 8192), performed 



on the most external points selected in the long trajectory of Figs. (11 and 



13). This graph reveals different structures near the plasma edge, allowing 



to choose the most external part of the trajectory, and a rational estimate oi 



the q- value of the external barrier : a KAM surface (or a robust Cantorus). . 34 



15 One of the last snapshot in the java animation. Bold points in this phase 



portrait (ip,9) represent the most recent iterated points along a single mag- 



netic line trajectory. Here, after filling the whole chaotic sea, the line is partly 



confined in the chaotic shell and then is sticking around the internal barrier 
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16 Farey tree indicating how the most noble number in a given interval between 
rationals (10/9 and 9/8, or 9/8 and 8/7) can be built with results N(l,8) or 
^(1,7)] ■ 

17 The upper Cantorus for L = 4.875/27T can be localized in the central region, 
between bold dots. Here we have plotted series of 2046 iterations along each 
of the island chains with rational q values corresponding to the following 
approximants to the upper Cantorus N(l, 7) : q — 17/15, 43/38 and 112/9S 
from above, and q = 9/8, 26/23 and 69/61 from below.| 

18 Part of the phase portrait of the standard map for K = Kg- Two series ol 

islands q = 8/5 and q = 5/3 can be seen, corresponding to two convergent^ 

towards the Golden number G. The last KAM with q = G is exhibited 

inbetween, and it appears as the last existing non-chaotized KAM surface; its 
initial conditions are I — 0.52969, 9 — 0.75. 

19 Calculated points (diamonds) of the perturbed g-profilc (around the magnetic 
axis) in the equatorial plane for L = 4.875/27T ~ 0.776, with the plasma ex- 

tending from X = —1 to +1, showing overall agreement with the unperturbed 

profile (dotted line), except for the the appearance of a local maximum on 

the magnetic axis (XMA).| 

20 Detail of the perturbed g-profile (measured around the magnetic axis) in the 
equatorial plane (circles) for L = 4.875/27T ~ 0.776, as compared with the un- 
perturbed profile (dotted line), showing the appearance of a local maximum 
on the magnetic axis, and of two minima. Horizontal levels represent (from 
top to botom) the g-values N(l, 7) = 1.131 of the upper Cantorus (which cor- 
responds approximately to the higher q- values of the bump), N(l, 8) = 1.116 
of the lower Cantorus, and N(l, 11) = 1.086 of the KAM barrier protecting 
the plasma core, along with the extension of these curves with the equatorial 
axis X. The value on the magnetic axis is quA — 1-1346 > iV(l,7) and lies 
at the center of a profile which is locally hollow at very small distances from 
the magnetic axis (not visible at this scale) 

21 The main noble transport barriers are represented in bold lines in normalized 

coordinates: " radius" = x = p/a = \fp and " angle" = 9. The robust KAM 

AT(1,11) separating the central protected plasma core (in white) from the 

chaotic shell. The two semi-permeable Cantori AT(1,8) and AT(1,7) form an 

internal barrier resulting in a very slow and intermittent motion towards the 

chaotic sea. The robust KAM torus on the plasma edge iV(4, 2) has been 

identified to have a vanishing flux 

22 Normalized histogram showing the occurrence H of sticking times longer than 
10 5 (in a run of 10 11 iterations with L — 4.875/27r), obtained by a sliding 
average, as function of the radial position ip of the visited island chain. The 
values of ip are taken within small intervals of 2.10~ 4 . This graph exhibits 
the frequent occurrence of long sticking times in and around the main rational 
chains : the two bumps correspond to widespread ip values of a line wandering 
in the main two chaotic zones: the chaotic shell (around ip ~ 0.1) and in the 
chaotic sea (around ip ~ 0.4 to 0.8) 

23 The figure visualizes the escape area from a bounded domain. The curve T(C) 
is the image of the closed curve (C) after one iteration by the map T . The 
domain Rc is the domain bounded by (C). The set T(C) n (C) consists in 
three unstable periodic points (namely U(l), U(2) and U(3)) and in three 
stable points (namely 5(1), 5(2) and 5(3)) 

24 Comparison of the decrease of the flux across convergent island chains towards 
the upper Cantorus q = N(l,7) (circles) and towards the KAM surface on 
the plasma edge q — N(4,2) (squares). The former is seen to converge to a 
finite value of the order of 10~ a , while the latter is observed with lower values 

to continue to decrease hopefully to zero (robust barrier). 



25 Average radius (a), mean square radial dispersion (b), average flux (c) and 




mean square flux dispersion (d) of an initial bunch of 500 lines at L = 4.875/2|7r 


followed during 10 y iterations 


26 Average radius (a), mean square radial dispersion (b), average flux (c) and 




mean square flux dispersion (d) of an initial circle of 5000 lines at L = 


4.875/2tt 




27 Average radius (a), mean square radial dispersion (b), average flux (c) and 




mean square flux dispersion (d) of an initial circle of 5000 lines at L = 5/2it 


28 Average radius (a), mean square radial dispersion (b), average flux (c) and 


mean square flux dispersion (d) of an initial circle of 5000 lines at L = 5.5/2tt 


29 Average radius (a), mean square radial dispersion (b), average flux (c) and 


mean square flux dispersion (d) of an initial bunch of 5000 lines at L — 6/2-71 




30 Average radius (a), mean square radial dispersion (b), average flux (c) and 




mean square flux dispersion (d) of an initial circle of 5000 lines at L = 9/2n 


31 Graphical presentation of the scaling dependence in the stochasticity param- 


eter L of the asymptotic coefficients of flux diffusion ~ L 2 (black circles) . 


of radius subdiffusion D x ~ L (grey squares), of average flux growth b ~ L 


(black triangles) and average radius growth a ~ L 1 ! 2 (black diamonds). For 


a lowest value Lq = 5.25/2-7T ~ 0.836 this graph represents these various co- 


efficients C in a logarithmic plot of log(C/Co) as function of \og(L / Lq) . The 


three straight lines correspond, from top to bottom, to the expected behaviors 


in L 2 , L and L Y I 2 , respectively 




32 Dispersion measurements of the radial position in the present situation of in- 




complete chaos are represented (lower black curve) for L — 5.5/2-7T ~ 0.876 (as 


in Fig.Dif 3b) with its asymptotic slope (dotted line) which exhibits an asymp- 


totic radial subdiffusion. In the domain between t — 10 3 and 10 4 (around the 


escape time t v , see Eq. 99), one clearly observe a superdiffusive regime. 


All these measured behaviors of MSD x (t) nevertheless appear at all rele- 


vant times to remain much smaller than the quasi-linear diffusion in complete 


chaos as represented by the upper straight line described by the classical eq. 


105 of Rcf.[69l (Rochester & Rosenbluth 1978) 





1 Introduction 



The ideal picture of perfect axisymmetric magnetic surfaces in a toroidal magnetic con- 
finement device like a tokamak, appears to be strongly modified either in presence of field 
inhomogeneities (e.g. divertors) or of tearing instabilities which result in the appearance of 
magnetic island chains, with a helical symmetry around the magnetic surfaces. The careful 
experimental analysis performed by N.J. Lopes Cardozo et al. [Q, with a very high spa- 
tial resolution on electron temperature T e and density n e radial profiles 0], up to a few 
times the width of an electron banana orbit, has shown that " small structures appear across 
the entire profile, with large magnetic islands occurring when the density disruption limit 
is approached" as a result of plasma filamentation He concludes that "the structures 
are interpreted as evidence that the magnetic topology in the tokamak discharges is not the 
paradigmatic nest of perfect flux surfaces, but more complex than that" . Similar results 
have been obtained on the TJ-II stellarator [Q where the very detailed structure of the T e 
profile measured along a chord has been observed with a A; -4 spectrum . 

In presence of several island chains, it is well-known that overlapping may occur J5|, 
resulting in the appearance of chaotic zones near hyperbolic points, and even chaotic seas 
with only island remnants. Equations describing magnetic lines in a torus are expressed in 
terms of the toroidal coordinate C along the line 0] . This variable C is usually interpreted 
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in analogy with a "time" variable, so that the Hamiltonian equations for magnetic lines 
are interpreted as "equations of motion". In a situation where island overlapping occurs, 
the classical picture describes a given magnetic line as "percolating" through the plasma, 
wandering in the chaotic sea, remaining almost trapped around island remnants ( stickiness ) 
and possibly reaching the plasma edge for increasing values of Q. In either case, this radial 
motion of the perturbed magnetic lines is responsible for an increased transport of particles 
and energy, since in the lowest approximation (with vanishing Larmor radii and vanishing 
magnetic drifts) charged particles just follow magnetic lines in the real time variable which is 
t ~ £ in this simple case. In the following description of magnetic lines, we mainly describe 
their so-called " motion" , for simplicity, in spite of the fact that a magnetic line of course 
remains static and is just followed in C, along the toroidal direction. 

It is known that a perturbed magnetic field configuration may sometimes be associated 
with the appearance of transport barriers (TB) characterized by a jump in the slope of 
density or temperature profiles, and a strong anomalous transport in the outside zone. 
TB could also appear in simulations due a strong velocity shear flow ("zonal flows"), even 
without magnetic turbulence ]tJ . In Ref. Q it is recalled that a strongly reduced ion thermal 
conductivity has been obtained in various tokamaks, which remains at the neoclassical level 
over part or even the entire plasma. This is due to the presence of an ion TB, associated 
with a strong velocity shear and a reduction of the density fluctuation level, which is not the 
case for electron TB observed due to the thin structures measured on T e |p[ . Such structures 
could be associated with alternative layers of good and bad confinement, localized between 
local barriers. It is important to stress the fact that transport barriers have been shown to 
be coupled to the safety factor profile (q-profile) and to exist also in ohmic plasmas. 

Experiments in Tore-Supra, JT-60U, JET, TFTR and other tokamaks (see review in Ref. 
||) clearly exhibit the influence of the safety factor profile (q-profile) on the appearance of 
internal transport barriers (ITB's) in tokamaks. For this reason we study here TB in a 
purely magnetic description. Electron ITB has been maintained during 2 s. in Tore-Supra 
Generally ITB's are obtained in presence of a reversed magnetic shear, i.e. with a 
q -profile which presents a local maximum near the magnetic axis, and a minimum q m i n 
typically at a normalized radius of 0.3 - 0.4, and then a regular increase towards the edge 
of the plasma. It is generally believed that ITB's appear around rational g-values, and may 
even follow the time-evolution of a given magnetic surface q = 2 for instance in JET lliof . 
Such ITB's may have a finite width and one finds an "ITB layer" [jlT) , defined as a thin 
layer with large T e gradients S inside the "ITB foot". 

On the other hand, even with a monotonia q -profile , a reduced heat diffusivity x nas 
been observed in the core region of the plasma, leading to the idea that ITB may appear 
even without reversed magnetic shear [fl2|| . This is the simple case we will mainly consider 
here. 

The experimental observation that ITB appear "near" rational surfaces may seem sur- 
prising from a theoretical point of view. Rational surfaces are not densely covered by mag- 
netic lines and thus are the most sensitive ones to plasma instabilities. It is a generic property 
that rational surfaces, on which field lines close back on themselves after a finite number 
of toroidal turns, are topologically unstable |jl3). Irrational surfaces, on the other hand, 
are covered by a single magnetic line, in an ergodic way, and appear to be more resistant. 
In the same way, the appearance of large scale chaotic motion is described in dynamical 
systems theory by successive destructions of KAM tori, and their transformation into per- 
meable dense sets, called Cantori. The most resistant KAM torus in a chaotic system is of 
course of crucial importance since it is the last inner barrier preventing large scale motion 
when the stochasticity parameter is increased. In simple cases it corresponds to a rotational 
transform (winding number) given by the "most irrational" number, the Golden number. 
From theoretical grounds, one can thus expect that irrational surfaces are more resistant to 
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chaos 14 and thus more likely to form a transport barrier, if any. In either case, stability 
or chaos, the breakup of magnetic surfaces "is a problem of number theory" fllq ], as will be 
verified here again. 

In the present work we are mainly concerned with the location of ITB's in a mono- 
tonic g-profile which appear naturally in a realistic model for toroidal magnetic lines, called 
"tokamap" ]l6| |. The unusual point we find here is that the magnetic perturbation, respon- 
sible for the appearance of the magnetic island chains, not only creates a non- vanishing 
Shafranov shift of the magnetic axis as expected, but also build a locally non-monotonic 
q-profile in the equatorial plane, with a spontaneous local maximum q max on the magnetic 
axis, and two local minima. The positions of the ITB found here appear to correspond rather 
exactly, in the perturbed q -profile , to "most noble" values of q which are the "next most" 
irrational numbers beyond the Golden number |Q [J These will appear to correspond 
to the q— position of the internal barriers of the tokamap. 

It has already been shown that the magnetic surface corresponding to a g-value given 
by the golden KAM torus is not the most robust barrier in the tokamap fl6|| . In other 
systems too, the golden mean is not found to be associated with the last KAM curve and 
the transition to global stochasticity [|l9), ^0|. Most noble values of q are however the 
location where we may expect, from KAM theory pl[ that the most robust tori are finally 
destroyed and changed into hardly permeable Cantori, which are thus good candidates to 
be identified with ITB's. That is what we will check to occur in the tokamap. This result 
does not fully agree with the generally admitted idea that ITB's in tokamaks would always 
be associated with rational q- values, but we have to note that the two Cantori forming the 
barrier found below are nevertheless observed on both side of a low order rational. 



Very interesting theoretical models based on transport across chaotic layers and internal 
barriers have been proposed to explain precise measurements performed on the RTP toka- 
mak. A model of radial transport in a series of chaotic layers has been developed |2^] where 
the standard magnetic equilibrium, with monotonously increasing q— profile, is perturbed 
by small closed current filaments: a number of filaments are localized on low order ratio- 
nal q— values, with suitably chosen values for their finite width and current. These current 
filaments break the topology of nested flux surfaces and are of course responsible for the ap- 
pearance of magnetic islands and chaotic regions. Test particle transport is computed and is 
found to be subdiffusive in such perturbed magnetic field, with a mean square displacement 
growing like the toroidal angle to the power ^ ]6^| . 

A model for inhomogeneous heat transport in a stratified plasma has also been devel- 
oped. Electron heat transport might of course be locally enhanced across each chain of 
magnetic islands (corresponding to rational q- value) and this could cause the appearance of 
the plateaux observed in the temperature profile M at rational q-values, explaining some 
jumps in the slope of the temperature profile. Such stochastic zones around low order ratio- 
nal chains, could also be limited by permeable Cantori. A simple analysis of heat transport 
in an inhomogeneous stratified medium shows that the measured values may deviate dra- 
matically from simple linear averages [|): the global transport description should take into 
account " insulating" regions but can ignore " turbulent" regions of high diffusivity. In such 
models for electron heat transport [^3| [Q, a number of transport barriers, with suitably 
chosen width and local heat conductivity, are assumed to be localized on surfaces with low 
order rational q~ values. As a result the large electron heat conductivity, assumed to be 
constant in the conductivity zone, presents a series of depletions {"q-comb model"). This 

1 \t is well known that the continuous fraction expansion of the golden number (the "most irrational" 
number) G = (y/5 + l)/2 = 1 + 1/ [1 + 1/(1 + ...)] = 1.618033989... is [1, 1,1, 1, 1....]. By changing the first 
1 to the left into an integer i > 1, one simply add unities to the golden number : G+ (i — 1) = [i, 1, 1, 1, 1....]. 
By changing the second 1 to the left into an integer n > 1, one obtains the "most noble" numbers N(i, n) = 
[i, n, 1, 1, 1....] which are the next most irrational numbers beyond the Golden one. 
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model succeeds to reproduce the changes observed in the temperature profile when scanning 
the deposition radius p dep of electron cyclotron heating (ECH) from central to far off-axis 
deposition. The main feature of these measurements p5| is the discontinuous response of 
the T e profile to a continuous variation of the deposition radius p dep : five plateaux, in which 
T e is rather insensitive to changes in p depl are separated by sudden transitions occurring for 
small changes of p dep . 

In all these models, the perturbations are assumed to be localized around low-order 
rational q— values, and even if several parameters have to be adjusted, the assumptions of 
these successful models are fully compatible with the result obtained here, i.e. the 
magnetic structure of an ITB as being composed of two noble Cantori. 

In Section (||) we present this simple Hamiltonian twist map ("tokaMAP" [jl6|) which 
has been proved to describe toroidal magnetic lines in a realistic way for tokamaks. Its 
derivation is summarized. We restrict ourselves mainly to the case of a monotonic q -profile 
in the unperturbed configuration. We derive some results concerning the localization of the 
fixed points, and determine the Shafranov shift. The bifurcations determining the number 
of fixed points are recalled, with an interesting example related to Kadomtsev's mechanism 
of sawtooth instability. 

Individual magnetic lines are calculated in Section (|J) for very large numbers of iterations, 
and a threshold region of the stochasticity parameter L is found, above which most lines from 
the central region actually reach the edge of the plasma {global internal chaos). The motion 
of a single magnetic line is found to be intermittent, with very long periods of trapping 
in different regions. The motion can indeed be localized between some layers separated by 
Cantori which correspond rather precisely to noble numbers in the profile of the perturbed 
q- values around the magnetic axis. The spatial localization of these barriers is discussed in 
Section ([|) . We also present in Section (||) the calculation of the flux through the Cantori 
barriers. A Cantorus is a fractal set of points, of fractal dimension zero in the poloidal plane 
|2^| (or of dimension 1 in the torus: a single magnetic line). Cantori are known to generally 
represent local permeable barriers in dynamical systems. 

In Section (|) we introduce an set of magnetic lines starting from a very small region 
(or a constant initial radius) and perform averages over this ensemble of lines. Iterations 
of the tokamap are interpreted in terms of "time evolution" describing the toroidal motion 
of a magnetic line. The average radial motion of the lines is described by calculating the 
"time"- dependent average radius and average poloidal flux 

r m (t)=<r(t)> , 1> m (t) =< j>(t) > (1) 

reached by the lines at each time, along with the mean square deviation of the flux coordinate 

-0 ~ r 2 , 

MSD4t)=<[i,{t)-4, m {t)] 2 > (2) 
and the dispersion of the radial coordinate with respect to this average radius 

MSD rm (t) =< [r(t) - r m (t)] 2 > (3) 
Here the brackets (...) indicate an average over the initial conditions. 

The time dependence of these three quantities are analyzed and the existence of an 
asymptotic regime is exhibited: in this regime we observe a diffusion of the flux coordinate 
(the variance grows linearly in time) , but a subdiffusion of the spatial motion (the variance 
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MSD rm (t) grows like t 1 ^ 2 ), and a still slower behavior of the average radius, which grows 
asymptotically like r m (t) ~ i 1 / 4 . Simple scaling laws are found to describe the dependence 
of the corresponding "diffusion" coefficients as function of the stochasticity parameter L. A 
quasilinear scaling in L 2 is obtained for the flux diffusion MSD^ , similarly to what is found 
in the standard map [^J, |2^|, for which additional oscillations deeply modify this simple 
quadratic growth. As expected, the subdiffusive radial motion MSD rm (t) has a coefficient 
characterized by a scaling in L 1 , while the slower radial motion r m (t) is characterized by a 
scaling in L 1 ^ 2 . 



We finally discuss the order of magnitude of the relaxation time, the time necessary to 
reach this asymptotic regime and discuss which regime could describe the magnetic line 
motion before reaching the plasma edge. The conclusions are summarized in Section (|7|) . 
Partial results have already been presented in a EPS-ICPP conference and published in p9f . 



2 Hamiltonian map for toroidal magnetic lines 

The unperturbed magnetic field realized in tokamaks is ideally represented by a set of nested 
toroidal magnetic surfaces wound around a circular magnetic axis. The condition V.B = 
allows us to express the magnetic field in the Clebsch form || in terms of the dimensionless 
toroidal flux ip and poloidal flux F. 

B = Vi/jxV9-VFxV( (4) 
2.1 Equations of motion for magnetic lines 

We use traditional toroidal coordinates (ip, 8, £) where 9 and £ are the poloidal and toroidal 
angles, respectively, and ip is the flux coordinate. From (Q) the "equations of motion" for 
the magnetic lines are easily derived: 

<hp_ _ 8F_ <W _ d£ 

d( ~ ~~80 ' ~aX, ~ dtp ^' 

These equations obviously have a Hamiltonian structure: the toroidal angle £ plays the role 
of " time" , and the poloidal flux F the role of the Hamiltonian. 

In the unperturbed case, F is simply a "surface function" F = F(ip) which represents 
an unperturbed Hamiltonian with one degree of freedom and corresponds to an integrable 
system: 

where 

Wty) = 1/<Z(V0 = tW/27T (7) 

is the winding number, the inverse of the safety factor q(ip)- Here the action variable ip ~ r 2 
labels the magnetic surfaces, it is canonically conjugated to the angle variable 8 . 

When a magnetic perturbation is applied, due to internal factors (instabilities, fluctua- 
tions) or to external causes (imperfection, divertor coils...), the poloidal flux F becomes in 
general a function of the three coordinates: 

F(iP,d,O=F o ^) + L6F^,0,() (8) 
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where L is the stochasticity parameter (). The field line equations become 

aQ 09 dQ dip 

corresponding to a Hamiltonian dynamical system with 1^ degrees of freedom, generically 
non integrable. The perturbation is responsible for the appearance of chaos. 

2.2 Hamiltonian maps, twist maps and standard map 

In order to avoid long symplectic integration in computing the magnetic line motion from 
(||) , discrete iterative maps have been introduced, specially to study the plasma edge. Many 
examples in the literature have been quoted in |L6|. An explicit iterative two-dimensional 
map consists in discrete transformations of the form 

1>u+i = p(1>M > ^+i = e(WA) (io) 

where v is a non-negative integer which represents physically the number of large turns 
around the torus, and where the functions P(^>„, 9 V ) and Q(ip v , 6 V ) are explicit in the "pre- 
vious" values ip v an d 9 V . 

Such transformations must of course conserve the Hamiltonian structure of the equations 
(|^) : the model should be a Hamiltonian map (i.e. area-preserving or symplectic) and 
therefore the transformation (|Io| ) has to be a canonical transformation of the canonical 
variables (tp, 9). In order to derive the map, one thus introduces a general generating 
function F(tp v+1 , 9 U ) which allows to write down the map as 

~ Ml ' ^+1 - — MT^l ( n ) 

which is in an implicit form. (Note that we choose here a generating function of the new 
momentum ip and the old angle 9, but the inverse choice is also possible and would lead to 
another family of maps). 

Other choices are possible for the same twist map, for instance with another kind of 
generating function F a (9 l ,,9 l ,+i) called the action generating function used in Section (JbJ), 
from which the map can be written as 

, _ 8f.(»„»,+i) a F„(fl,,,fl„ + i) fl2 x 

The relation between F and F a is thus (from (|ll]) and (|l2|)): 

F a (9v, 9 v+ i) = ip u+1 (9 v , 9 u+ i).9 v+ i ~ F \%p v+l (9 Vl 9 v+ i), 0„] (13) 

(seeRef. @). 

2.2.1 General form of a twist map 

In order to proceed, the following general form of the generating function has been chosen: 

F(ip v+1 ,6 u ) = i> v+l .e v + FoWv+i) + L 5F(ip v+1 ,9 u ) (14) 

where L is the perturbation parameter and i'bWv+i) the unperturbed term taking into 
account the winding number (see(Q)) 

dtp q{ip) 
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As a consequence of ( pi] ) the map takes the following form 

VVfX = i> v + L h(l/) u+1 ,8 u ) 

e u+1 =e v + w(xp v+1 ) + l jtyu+M (16) 

which is an Hamiltonian form because the two functions h and j defined by 

h(^ 1/+1 ,6 v ) = , 3(.i> v+1 ,0 v ) = (17) 

automatically insure that 

d h(i/) v+1 ,9 u ) d j(^ v+ i,O u ) 



(18) 



Equations (|16| ) has the form of a general Hamiltonian map, from which simple cases can be 
recovered. 



When variables are clearly separated in (|16|) i.e. for h = h(9„) and j = one 



finds a general twist map in the case where the profile q(ip) is monotonous (see section [5.1 
below) . 

2.2.2 The standard map 

The Chirikov- Taylor standard map J27[ @] belongs to that family and corresponds to 

h{9) = -sin 2it9 , j(i/>)=0 , W(4>) = ip (19) 
and has the following form 

= i> v - Ls\d-2tt9 v , 0„+i = 9 V + ip v+1 (20) 

corresponding to 

F (^) = iV 2 , W(lM)=-£coB(27r0) (21) 
The generating functions for the standard map are thus 

F(fl> v+1 ,9 v ) = i> v+1 .9 v + - |- cos(27T0„) (22) 

and 

F a [9 vi 9 v+ i) = \\B„- 0„+i] 2 + cos(2tt^) (23) 

2 Z7T 

2.2.3 The Wobig map 

The Wobig map ]3l| on the other hand corresponds to : 

h(tp,6) = -<0sin27T0 , j(ip) = cos27r6> , W(i/>) = tp (24) 

and has the form 

=ip u - Lipv+i sin27rfl„ , 0„+i = + Vv+i - T^rr cos27r^ (25) 
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corresponding to 

F {ib) = \ib 2 , 8F{i>,$) = -^coB{2ir0) 
The generating functions for the Wobig map are thus 

1 2 L 

F(tb v+1 ,9 v ) = ip v+1 .0 v + 2^+1 - 2tt^ 1/+1 cos ( 27r6 ' l/ ) 



and 



F a (9 v , 6 v +i) — 



-i 2 



/+1 - #y + 7T" COS(27T0„) 
Z7T 



(26) 



(27) 



(28) 



These last two maps are actually not suitable to represent magnetic lines in a tokamak 
first because they do not insure that a nonnegative value of ib ~ r 2 remains nonnegative after 
iteration (as it should to represent a real value of the radial position) and, second, because 
they do not involve any realistic profile of the winding number W(ib) (they correspond to a q 
-profile everywhere decreasing: q ~ 1/r 2 ). In order to satisfy these two necessary properties, 
another model is described in the next section. 



2.3 The tokamap 



A specific model, the "tokamap" has been derived in fllq | from general properties of Hamil- 
tonian twist maps. The advantage consists in succeeding to describe the whole body of 
a tokamak plasma, including chains of magnetic islands in a realistic way. This map de- 
scribes the basic motion of the magnetic lines in the two dimensional poloidal plane (ib,0) 
by the winding number W(tp) which is here modified by an additional contribution from the 
magnetic perturbation. 



The general expression of the tokamap results from the following choice in the generating 
function (14) : 



1 ib 

5F(ip +1 ,0 U ) = — — cos2tt9 v (tokamap) 

2tt 1 + ip v+1 



(29) 



which involves an additional ip dependence as compared to Eq.(Pq). We thus consider the 
following generating function for the Tokamap (see (|l4|) and (p9|)): 

ib 

F(tp u+1 ,6 u ) = ib u+l .6 v +FqU> v +i) ~ L - vJ \ x cos27r6> ;y (tokamap) (30) 

1 + Wv+l 

This immediately leads to the following implicit form of the tokamap (use (flTl)): 

1> v+1 =ib v - L ^ Bm(2n0 v ) (31) 
1 + Wu+l 



e v+ i = e v + W{^ u+ {)-——— ^cos(2^) (32) 

2tt (1 + Vv+i) 2 

where 9 denotes the poloidal angle divided by 2-7T . In this nonlinear map a unique root is 
chosen for Eq. fl3l|) : 

1>u+i = \ (P(^, 0u) + vW„A)] 2 +4</v) (33) 
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where 



P(i>,6) = ip-1-Lsm(2n9) (34) 

This explicit map ([3l]-|34|) has been shown to be compatible with minimal toroidal geometry 
requirements. The polar axis {tp — 0) represents the magnetic axis in the unperturbed 
configuration. This map depends on one parameter, the stochasticity parameter L (R), and 
on one arbitrary function W(tp) which is chosen according to the q -profile we want to 
represent. For increasing values of L, chaotic regions appear mostly near the edge of the 
plasma. 

It is a simple matter to check that the symmetries of the tokamap imply that, for negative 
values of L, the phase portrait is the same as for positive values, but with the simple poloidal 
translation : 9 => 9 + 1/2 which means that the original phase portrait for 9 = to 1 is 
recovered identically but for 9 = —1/2 to +1/2 . 

It has been shown that, contrary to what occurs in the Chirikov- Taylor standard map 
pTj or even in the Wobig map fl3~i| , a nonnegative value of tp always yields a nonnegative 
iterated value, insuring that the radius remains a real number (see Eq.( |123| )). In the domain 
L < 1, however, the Wobig map also conserves nonnegativity of ip . 

The tokamap has recently been deduced in a very different way ]32| , as a particular case 
of a particle map describing guiding centre toroidal trajectories in a perturbed magnetic 
field given in general by (Q) , and to lowest order by the standard magnetic field model J33| , 
l3l: 



B 



1 + Et x COS ( 



X 

q(x) _ 



(35) 



where Et = a/Ro is the inverse aspect ratio of the torus, i?o the large radius of the torus, 
Bo the magnetic field strength on the axis and x ~ p/a the normalized radial coordinate of 
the plasma of small radius a . This magnetic field ( J35[ ) is checked to be divergenceless. By 
comparing the poloidal flux expressions in the general Clebsch formula (|J) and in the above 
standard model ([35]), it is easy to deduce that the toroidal flux tp in this case is given by 

^ = y (36) 
in the dimensionless units used in the tokamap. 

In [ p2| canonical coordinates for guiding centre have been derived, allowing for the sym- 
plectic integration of the equations of motion and the derivation of a Hamiltonian map 
for guiding centre in a perturbed toroidal geometry. As a particular case, the tokamap is 
deduced from this particle map when one applies a simple m = nonresonant magnetic 
perturbation. In order to describe magnetic line motion only, the magnetic moment is con- 
sidered to be zero, and one keeps terms to the lowest order in the inverse aspect ratio £t- 
Moreover, in order to go from the time dependence of the particle trajectory to the toroidal 
(■-dependence of the position (ip,9) of the magnetic line, all equations are simply divided 
by the equation for d(/dt . As a result, Eqs.(|l]-|34|) are exactly recovered. This deriva- 
tion also allows to write down explicitly the form of the magnetic perturbation involved in 
the tokamap : The generating function ( |l4] ) takes into account a magnetic divergence-free 
perturbation <5B with the following components: 



8 J?± = T X - eT gjgMlf) f o ?1 

B ~h(x,9)' dij [ ' 



2 The stochasticity parameter used here for convenience is related to the parameter K of the original 
reference by Balescu, Vlad & Spineanu [16] by: L = K/2n. 
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and 

B 'x.h{x,ey de [ ' 

where h(x, 9) = (1 + et-x. cos 9) represents the toroidal effect. 

A different but related method of derivation of symplectic maps has been used by Ab- 
dullaev §|. 

2.3.1 Safety factor profile 

In order to take magnetic shear into account, a specific q -profile has to be introduced to 



describe the unperturbed equilibrium. We consider here the same profile as in Rcf. |16 



S#) = — ; — — o- (39) 

which corresponds to a classical cylindrical equilibrium [ ^6| in which the value on the axis 
(if) = 0) is given by the parameter w : 

<z(V> = 0) = - (40) 



w 



and on the edge (ip = 1) by : 



q{i> = 1) = 4g(0) (41) 



which is four times the central value. The derivation of this q -profile ( J39| ) is presented in 
Appendix A. 

2.3.2 One example with reversed magnetic shear : the "Rev-tokamap" 

If a non monotonic q -profile is used instead, a reversed magnetic shear is introduced in 
the unperturbed magnetic field, and it has been shown that the above mapping becomes a 
nontwist map : this "Rev-tokamap" |S7| has very different properties namely because the 
Kolmogorov-Arnold-Moser (KAM) theorem does not apply anymore. It has been found |37j 
that a critical surface appears in the plasma near the minimum of the q -profile , separating 
an external, globally stochastic region from a central, robust nonstochastic core region. Such 
a phenomena of " semiglobal chaos" has been shown to be analogous to the appearance of 
ITB in reversed shear experiments. Later, similar theoretical results have been obtained in 
the Lausanne group ]3§[ ] in a reversed shear TEXTOR equilibrium; by studying a symplectic 
perturbed map for the magnetic topology, they confirmed that the transport barrier is indeed 
localized near the minimum of the q -profile . 

Coming back on the reversed shear model of Ref. [B7| we will show below in Section 



(4.5.2) that the ITB found in that work is also localized with a good precision on a g-value 



which will appear to be a noble number. 



2.3.3 Fixed points and bifurcations 



The general structure of the phase portrait of the tokamap has been described in detail in 
p^ j. Fixed points of the map should not be confused with secondary magnetic axes inside 
the islands: the latter are not fixed but wander upon iteration from one island to another 
island of the chain and appear as periodic points. 
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Fixed points, defined by i/v+i = an( i 9v+i — ®v + n where n is an integer, can 
easily be calculated by using first Eq. ( pl| ) which yields either -0 = (the polar axis) or 
sin(27r0i/) = 0. This proves that all fixed points are either 

- on the polar axis (ip = 0), for instance 9 = 0.25 and 0.75 which are two hyperbolic 
points on the polar axis, and exist as long as L < 1, 

- or in the equatorial plane (9 = or 0.5). By using next Eq. (|3^), one finds that the 
values of ip at the fixed points with 9 = are solutions of 

Fiftft, K, w, n) = Wty) - ^ + * - n = (42) 

while values of ip a t the fixed points with 9 = 0.5 are solutions of 

JM, K, w, n) = W(il>) + ^ + X - n = (43) 

It has been shown |p^| that the origin, the polar axis (ip = 0) is a fixed point as long 
as L < 1. Fixed points have been discussed in great detail in including consideration 
of the "ghost space" ip < which appears necessary in order to check that the Poincare- 
Birkhoff theorem and the conservation of stability index (the difference between the number 
of elliptic and hyperbolic points which remains constant as w varies) are indeed satisfied. 
Let us simply recall here that bifurcations occur according to the value of w representing 
the central winding number (see ^0|) as compared to the following L-dependent parameters 

w m (L) = 1 — L (44) 

and 



w M (L) =2-w m = 1 + L >w m (L) (45) 

As long as w < w m = 1 — L, the only invariant point is the polar axis ip = 0, which obviously 
has to be interpreted in this case as the magnetic axis of the tokamak (see Fig. 15 in Jj"6||). 

For higher values of w, when w m < w < wm = 2 — w m , i.e. 1 — L < w < 1 + L, a first 
bifurcation has occurred, with the appearance of an elliptic point off the polar axis. This 
latter elliptic point has to be interpreted as the magnetic axis, displaced by the Shafranov 
shift |3£| , |Q . In the case w = 1 (to which we will restrict ourselves in the next Sections) , 
the position ip MA of this magnetic axis is obtained by using Eq. ( |43| ) as function of L as: 

F 2 (iP MA ,L,w = l,n = l)=0 (46) 

The curve giving ip = iPma{ k = 27rL ) is plotted in Fig.(@). 



In the absence of any perturbation (L = 0), the magnetic surfaces are centered around 
the polar axis. 

A first typical phase portrait is given in Fig. 5 of Ref. |ll| for the simple case where 
q(0) = 1, and L = 4.5/27T ~ 0.716 , hence w = 1 and w m < w. One can see the formation 
of a chaotic belt which surrounds several large islands chains. This belt is clearly confined 
between two surfaces which were interpreted as KAM surfaces in this case. Several other 
examples are given below (see FigsJ^-||). We will discuss why, for slightly higher values of i, 
such boundary surfaces can be identified as fractal barriers or Cantori, which can be crossed 
by a magnetic line after a very long time. 
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Figure 1: Analytic solution for the Shafranov shift: position ip of the magnetic axis as 
function of K = 2irL in the case w = 1. For L = 4.875/27T this formula yields ip(L = 
4.875/27r) = 0.07492 . Here r — \p§ represents the radial coordinate, with r = y/2 on the 
edge. 

2.3.4 Second bifurcation for q(0) < 1 

A second bifurcation occurs for still higher values of w when w m < wm < w, i.e. when 
1 + L < w which implies q(0) < 1 : the magnetic axis remains displaced from the origin 
and corresponds to the center of a m — 1 island in the equatorial plane at 9 = 0.5, resulting 
in the appearance of an hyperbolic point ("A" point") for 9 = with a separatrix enclosing 
the origin. The latter appears to be the "main" magnetic axis, "expelled" by the presence 
of the m = 1 island. This can be seen in Fig.(Q) which presents the tokamap phase portrait 
for L = 1/47T, w = 1.05 (thus q(0) < 1 and w m < w) over N = 8193 iterations. 

This example shows the final situation after a big island (m = 1) has grown and expelled 
a small region around the original magnetic axis near the polar axis; this occurs when the 
unperturbed q- value on the axis is smaller than unity, q(0) = 1/w < 1, allowing for the 
existence of a surface q = 1 inside the plasma. It has been stressed that this process occurs 
as a result of a reconnection of magnetic lines, and is typical of Kadomtsev's theory of 
sawtooth instabilities. The position of the m — 1 magnetic axis (" O point" ) , as well as the 
position of the hyperbolic X point can be calculated analytically from Eqs. ( |42"[ in 
agreement with Fig. 17 of Ref. Jl(|. (We note that two misprints appeared in the text 
of that paper for the exact numerical values which should be written ipl = 0.040 of the O 
point, and ip2 = 0.024 of the X point, as indicated by the Fig. 17, instead of 0.179 and 0.167 
written in the main text). 

These points can be calculated as follows. For wm < w the position tp H of the hyperbolic 
X point is given by the solution of Fi(ip H , L, w, 1) = (with 9 — 0) as defined in Eq. ((42|). 
For a given w, this solution exists only for the smaller values of L. For instance for w = 1.05 
it exists only for K = 2-kL < 1.97392 , which means that for q(0) < 1 this X point only 
appears for small values of the perturbation parameter. 
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Figure 2: Tokamap phase portrait for L = 0.5/2n with q(0) < 1 (w = 1.05 and N = 8190). 
The polar axis (0, 0) is the main magnetic axis of the plasma core. The position of the 
elliptic point in the m = 1 island is if = 0.040246, 9 — 0.5 and the hyperbolic point is 
located in ip = 0.024481, = 0. 
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The position ip E of the m = 1 magnetic axis (elliptic or O point) is given by the solution of 
Fzi^Ei Li w, 1) = (with 9 = 0.5). We remark that this last solution ij) E {L) does not vanish 
for L =>• 0, which would seem to indicate the existence of a fixed point [ip E (L = 0), 6 = 0.5] 
out of the polar axis '0 = 0, even in the unperturbed case where all magnetic surfaces are 
circles centered around the origin... In this case we have wm(L = 0) = 1 < w which means 
q(0) < 1 and this fixed point [ip E (L = Q), = 0.5] is simply a standard point of the circular 
surface q — 1 : this is just a special feature of this unperturbed surface on which every point 
appears as stationary since it exactly comes back after one toroidal rotation (one iteration) . 
This specific point with 9 = 0.5 in the unperturbed case is however important since it acts 
as a seed for the m = 1 island which appears when L becomes different from zero. We note 
that this position ip E (L) is such that q(i/j E (L = 0)) = 1. 

In order to simplify the discussion about dispersive motion of the magnetic lines in the 
tokamap, we will restrict ourselves in the following Sections to the case 

w = 1 (47) 

which means q(Q) = 1 . The unperturbed q -profile is a continuous and monotonous function, 
with a value growing from 1 in the center to 4 on the edge of the plasma, a rather standard 
profile in most ohmic discharges. In this case for any L > 0, we have from Eq. ( p4| ) 
w m (L) = 1 — L < w, i.e. we are beyond the first bifurcation and the position of the 
magnetic axis is given by Fig. ([[]). 



3 Individual trajectories: tokamap phase space portrait 

In this section we first present the phase portrait of the tokamap in the case w = 1, q(0) = 1 
for increasing values of the stochasticity parameter and then study very long trajectories. 



3.1 Increasing stochasticity 

For increasing values of the stochasticity parameter the tokamap exhibits all the main fea- 
tures of chaotic systems, in a way which is very realistic for tokamaks. 



3.1.1 Weak stochasticity regime : confinement by KAM tori 

For small values of the stochasticity parameter L = K/2ir << 1, most of the KAM surfaces 
are preserved and the phase portrait of the tokamap appears to be described by embedded 
tori, around a magnetic axis displaced from the origin (ip = 0, 9 = 0) by the Shafranov shift 
(46), as seen in Fig.([|). 



In this case all magnetic lines remain confined, no magnetic island is seen. From the 
position of the shifted magnetic axis, we note that, for positive values of L, the weak field 
side of the torus is in the direction 9 = 0.5 . 



3.1.2 Appearance of island chains on rational g-values 

For an increased value of L, regular chains of magnetic islands appear; their number m and 
the order in which these islands are visited by a magnetic line allow us to deduce their q- value 
in a simple way. For instance a chain of m islands visited one by one at each iteration, in the 
direction of increasing values of 9, has a rotational transform or winding number t(-0)/27r 
equal to 1/m., and a q- value q = m. These are the main island chains, as seen in Fig.(|4]). 



A direct measurement of the exact value of q can be performed, specially for almost 
undestroyed magnetic surfaces, by computing the average winding number along a trajectory, 
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Figure 3: Regular KAM magnetic surfaces with a Shafranov shift for a small value of the 
stochasticity parameter L = 0.08, with w = 1. The angle 9 is represented modulo 2n. 
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Figure 4: Regular KAM tori for L = 0.4, w — 1, and main rational island chains corre- 
sponding to q = 5/4, 4/3, 3/2, 2, 3 and 4 
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see (0). For small circles around the magnetic axis (which do not encircle the polar axis), 
a change of coordinate is necessary to evaluate correctly the rotation around the magnetic 
axis. On the other hand, for all surfaces encircling the polar axis, the knowledge of the 
iterated values 9 T allows us to calculate the following surface quantity (average on a given 
magnetic surface) 



which is the average increase of the poloidal angle over a large number N of iterations. The 
q— value of the corresponding trajectory is immediately obtained from (|t]). 

3.1.3 Overlapping chains, secondary islands and chaotic regions 

For an increased stochasticity parameter L — 4.875/27T ~ 0.776 , the same initial conditions 
as in Fig.(|j) actually describe a wide chaotic zone surrounding a central part with regular 
KAM tori : see Fig.(|). 

The main rational chains q — 3/2, 2, 3 and 4 are partly stochastizcd. 

We show on Fig.(^) a small part of the phase portrait obtained by following 6 trajectories 
along 16 380 iterations: the main chain of primary islands is observed, with intermediate 
higher rational chains, surrounded by stochastic regions in the vicinity of the hyperbolic 
points (as usual in nonlinear dynamical systems). 



3.1.4 Island remnants in the chaotic sea 

The general aspect of the phase portrait is determined mainly by the q -profile . The 
structure of the stochastic sea is analyzed in terms of the various island chains. The islands 
are partly "destroyed" by chaotic regions on their edge, so that we mainly observe "island 
remnants" which form a "skeleton" of the phase space. Several island chains among the 
largest, which belong to the "dominant" classes n = m and n = m — 1 with 



are represented in Fig.(Q). This skeleton structure is then invaded by the chaotic sea (not 
represented here). 

This last series n = m — 1 has some importance here : it corresponds to chains of m 
islands visited one by one at each iteration, but with decreasing values of 9, thus with a 
rotational transform l/2tt equal to — 1/m (modulo 1), or equivalently to (m — l)/m, hence 
q = m/(m — 1). 

This global island structure (or geographic chart of the stochastic sea or "skeleton") is 
only slightly perturbed for larger values of L. In Fig.(^) we have represented for L = 
5.25/27T ~ 0.836 a series of small island remnants appearing as resistant KAM surfaces 
around the secondary magnetic axis inside the islands ("vibrational KAM" fl4|| ). From top 
to bottom we observe 

- chains of island corresponding to [m = 4, q = 4], [m = 3, q = 3], [m = 2, q = 2], 

- along with smaller islands: [m = 3, q = 3/2], [m = 4, q = 4/3], [m = 8, q = 8/7], 

- a stochastized region around the [m = 9, q = 9/8] island remnants, 

- and also [m = Yl,q = 12/11 = 1.09]. 

All these chains are seen to encircle both the magnetic axis and also the polar axis (since 
they cover the whole interval of 9). We also represented two chains encircling the magnetic 
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Figure 5: With the same initial conditions as in Fig. 4 (w = 1) with a higher value of 
the stochasticity parameter L — 4.875/27T ~ 0.776, regular KAM tori are observed in the 
central part, while the main rational island chains corresponding to q = 3/2, 2, 3 and 4 are 
stochastized. 
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Figure 6: Partial vue of the phase portrait with six trajectories for L — 4.875/27T ~ 0.776, 
which shows the structures of rational chains of island remnants and the stochastic web 
expanding from the hyperbolic points. 
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Figure 7: Part of the skeleton of island remnants in the phase portait (r = \/%p = as 
function of 8) for L = 4.875/27T ~ 0.776 : rational chains can be seen corresponding to 
q = 4, 3, 2, 3/2, 4/3, 8/7, 9/8, 12/11 . 
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Figure 8: Several trajectories for L = 5.5 /2tt ~ 0.875.. , yielding a central protected plasma 
core, surrounded by chains of regular islands and KAM surfaces, and, outside, a stochastized 
zone allowing the line to escape from this central chaotic shell. 

axis only : [to = 13, q= 13/12 = 1.083] and [m = 11, g = 11/10 = 1.1] which is inside the 
former ones. 

From the respective positions of the two latter chains, we note the unexpected fact 
that, locally, the perturbed q-value is actually growing towards the magnetic axis, a very 
important feature for transport properties. This proves that a non-monotonous q -profile 
around the magnetic axis, and a negative magnetic shear has been spontaneously created by 
the magnetic perturbation. The precise perturbed q -profile is presented below in Section 
|J (sec Figsg 

3.1.5 Good confinement of the plasma core 

Between these island remnants are large chaotic regions. Of particular interest is a circular 
chaotic shell (or belt) around the magnetic axis. For instance with L — 5.5 /2tt ~ 0.875 
we have represented in Fig.(|J) several trajectories filling a circular shell surrounding the 
magnetic axis, and wandering in a chaotic layer around island remnants. We stress the fact 
that this chaotic central shell surrounds a regular central part (around the magnetic axis), 
which represents a quiet central plasma core protected from chaos, thus a good confinement 
zone. 



We note that the chaotic zone represented here has an obvious internal separation, with a 
visible trajectory passing very near r ~ 0, 9 = 0.25, where there is actually a fixed (unstable) 
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Figure 9: Phase portrait of the plasma edge, after complete destruction of the KAM surfaces. 
Several trajectories are chosen to represent the details of the four island remnants q = m = 4 
at a value L = 6/2n ~ 0.954... beyond the escape threshold. The plasma edge has become 
permeable and is strongly deformed as compared to the unperturbed circle '0 = 1. One 
remarks the presence of four satellites or daughter islands around each main island. 

hyperbolic point ]l6| ]. The separation is actually nothing else than the chaotized separatrix 
in the (ip, 9) plane, joining the two hyperbolic points located in rp = 0, 9 = 0.25 and 0.75. 
In a polar representation this separatrix is the circular magnetic "surface" encircling the 
magnetic axis and tangent to the polar axis rp = 0. 



3.1.6 Stochasticity threshold for escaping lines out of the plasma bulk 

The aim of the present work consists, first, to determine, for a given q -profile , the threshold 
domain of the values of the stochasticity parameter L for which the plasma boundary be- 
comes permeable, allowing the magnetic line to escape across this broken barrier to the edge 
and corresponding to a disruption of the plasma towards the wall. This happens between 



L = 4.875/27T ~ 0.776 and L — 5/2ir ~ 0.796. The next question is (see Sections 4.1 and 



4.3) : for lower values of L, in the confinement domain, which is the most robust barrier 
able to inhibit the motion of the magnetic lines up to the edge ? In other words : which is 
the most robust KAM surface inside the plasma, the last one to be broken ? Which is its q 
value ? 



For L = 6/2-7T we remark that the edge of the plasma has become permeable and is 
strongly deformed as compared to the unperturbed circle tp = 1. The four islands q — m — 4 
can be seen on Fig.(0), along with their four satellites or "daughter islands". This value of 
the stochasticity parameter is obviously larger than that of an escape threshold. 
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Figure 10: Single trajectory at L — 5.5/27T ~ 0.875 followed during N = 16380 iterations. 
After a long trapping period in the chaotic sea, this trajectory finally escapes out of the 
plasma edge. The initial conditions are r = y/Zip = 0.4 and 8 = 0.025.. 



By performing a very long iteration up to 2.10 9 time steps on a Alpha workstation, with 
a value L — 4.875/27T ~ 0.776, we did not reach the time where the particle could possibly 
escape from the plasma : up to this time the trajectory remains confined by what can be 
considered as a KAM surface. 

On the other hand, for slightly larger values of L of the order of 5/27T ~ 0.796, most of 
the magnetic lines are found to rapidly escape from the plasma (see for instance Fig.([To|) 
for L = 5.5/27T ~ 0.875). A threshold region of the stochasticity parameter has thus been 
found slightly below L — 5/2tt ~ 0.796. For larger L values, magnetic lines escape across 
the plasma edge even when starting from the central chaotic shell. 

3.2 Intermittent motion inside the confined plasma: crossing inter- 
nal barriers 

For values of L smaller than this escape threshold, magnetic lines are wandering inside 
the plasma. We mainly consider here the case L = 4.875/27T ~ 0.776 at which the edge 
barrier is not yet broken. We follow one single magnetic line along a large number of 
iterations, starting near the central region. With this value of the stochasticity parameter, 
the corresponding phase portrait is particularly rich, even with only one point represented 
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Figure 11: Phase portrait for L — 4.875/27T ~ 0.776 as drawn in coordinates 9 and r = y^tp 
by a single trajectory followed during 2.10 9 iterations, with one point represented every 10009 
iterations. The inner shell appears slightly darker, surrounding the plasma core (white) and 
the magnetic axis. Long sticking occurs here around the five small q = 5/2 islands (dark 
lines in the upper part of the chaotic sea). 

every 10009 iterations 0. 

The chaotic sea is represented on the Fig.(|ll]) where we recognize the protected zones 
corresponding to the island remnants with q = 4, 3, 2 and 1 (the central core) beside all 
other main rational chains. 

For this trajectory the sticking stage is particularly long around the boundary of the 
q = 5/2 island remnant, and lasts for several 10 8 iterations. Starting from points in the 
neighborhood, and performing a small number of iterations (N = 16380) we find that this 
dense stochastic zone takes a figure-eight form (which indicates that a period doubling has 
occurred) with 9 daughter islands around (see Fig.(|l2"|)). Each of the central elliptic points 
in these q = 5/2 islands has already bifurcated, giving rise to an inverse hyperbolic point 
and to two new elliptic points. 

It is most interesting to determine the time behavior of the magnetic line position in 
its complicated motion across the chaotic sea of the poloidal plane. The full information 

3 In selecting a finite number of points on the graph, in order to avoid a fully black drawing, we have 
represented in this case only one point every 10009, a large prime number to avoid lower order graphical 
stroboscopic effects (by using any multiple of 5 for instance we would have drawn only one of the five islands 
9 = 5/2 since the same magnetic line visits the neighborhood of these five islands one after the other - two 
by two - and comes back in the first island neighborhood after 5 iterations. In order to represent the whole 
set of 5 daughter islands around the two islands q = 2 we need to avoid multiples of 10, a.s.o... so that we 
choose a large prime number). 
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Figure 12: Detail of a small structure of the phase portrait around one of the five q = 
5/2 islands, where a long sticking stage occurs during 10 8 iterations along the trajectory 
represented in Figs. 11 and 12. 
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has been represented on an animated movie [ pd| . We have followed a very long trajectory 
on 2.10 9 time steps starting from ip = 0.0136125, 8 = 0.033. As discussed in Appendix 
A (Eq. ( |123| )), this initial condition corresponds to a radius r — 0.165 , slightly out of 
the central protected core of the plasma (we recall that in these notations, the edge of the 
plasma ip — 1 is at r = y2 ). From such initial conditions, the road is completely open up 
to the edge of the plasma, which is a resistant KAM barrier for this L value, but the path is 
far from trivial: the magnetic line describes a path "percolating" among the rich variety of 
island remnants which are known to form a hierarchical fractal structure. In that movie one 
remarks that the line remains for a very long "time" in a chaotic layer, between the central 
protected core and some visible transport barrier "around" q = 9/8, and then crosses this 
barrier to explore the upper chaotic sea (extending up to the plasma edge). Such transitions 
across the barrier occur repeatedly, inward and outward, in an intermittent way and after 
random residence times. Shorter periods of sticking inside this barrier are also observed. 

This behavior can also be represented in a (r, t) graph shown in Fig. ( |l3| ) where we present 
the time variation of the radial position : this behavior is not only stochastic, but presents 
different stages : 

- a first period of pseudo-trapping or temporary confinement in an inner shell below 
r^0.6, 

- an escape in an external shell, and a wandering stage in a wide chaotic sea, avoiding 
island remnants, up to r ^ 1.6, 

- a period of deep sticking around the five island remnants q = 5/2 (only four bands can 
be seen since two of the five islands overlap in radius), 

- a new wandering in the chaotic sea of the external shell, then a second stage of pseudo- 
trapping in the inner shell, 

- and so one and so forth... 

In other words the oscillating radial position of a single line is seen to wander randomly 
from the (inner) stochastic shell, to the stochastic sea (sometimes around q = 5/2 island 
remnants) with small downward peaks indicating excursion towards and inside the transport 
barrier, like flood tide and ebb tide. 

Such an intermittent behavior is typical of the phenomena of alternative trapping into 
different "basins" of a Continuous Time Random Walk (CTRW) Q, as previously observed 
|43[| in a stochastic layer of the Chirikov- Taylor standard map. One may wonder why these 
average residence times (trapping times) are so long in the present case. This could be 
related with the fact that the value L — 4.875/27T considered in Fig.(|l^) is not far from the 
escape threshold L ~ 5/2-7T, and that the latter could be the analogous of a critical point 
of percolation, in the vicinity of which cluster lengths and diffusion characteristic times are 
generally diverging as some inverse power of the deviation from the critical parameter value, 
with a critical exponent. 

This conjecture however remains to be proved. In the present problem, such long char- 
acteristic times seem to exclude the possibility to compute the histogram (the probability 
distribution of these residence times), which would need awfully long iterations times in 
order to obtain a good statistics. 

4 Localization of transport barriers 

The intermittent motion described in Fig. ([13]) clearly exhibits intermittent periods of con- 
fined motion between structures playing the role of internal transport barriers (ITB), the 
most resistant curves inside the confined plasma. It is interesting to identify such barriers 
and to note their positions along the q -profile . 
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Figure 13: Radial position of a single trajectory followed during 2.10 9 iterations, as func- 
tion of time for L = 4.875/2tt ~ 0.776,with initial condition tp = 0.0136125, 6 = 0.033 : 
different stages of trapping in different zones indicate an intermittent behaviour. Note that 
trapping stages are very long at may last for at least several 10 iterations for this value of 
the stochasticity parameter, near the escape threshold. (Such long times near the critical 
threshold for large scale chaotic motion could be reminiscent of critical phenomena. 
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4.1 Rough localization 



4.1.1 Central core barrier 

An initial period of pseudo trapping is observed, during which the trajectory remains in the 
inner shell of Fig. (pi]) for several 10 6 iterations. This shell is separated from the protected 
central core of the plasma by a very strong barrier on the plasma edge which has not been 
crossed during 2.10 9 iterations. This barrier can thus be considered as a KAM torus (or 
a very robust Cantorus if it could have been crossed by performing still more iterations). 
By analyzing the innermost points reached by the long trajectory, and performing short 
iteration series on several of them, we can determine which part of the trajectory is a good 
candidate to localize the barrier. A rational estimate of the value of the safety factor can be 
determined on this trajectory by carefully analyzing the (almost) periodic repetition of the 
variation of the poloidal angle 9 in time, along with an estimation of the average poloidal 
rotation at each iteration, or by a direct calculation of the winding number ( fl8| ) after a 
possible change of coordinates for innermost trajectories not encircling the origin. We find 
that a very rough estimate of the inner barrier protecting the plasma core from invasion 
from the inner shell is characterized by: 
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qci = = = 1-080 (50) 



4.1.2 Internal barrier preventing outward motion in the chaotic shell (lower 
Cantorus) 

In this inner shell, the outward motion is limited by a series of points which has been 
analyzed in a similar way. Measurement of the g-value of the outermost trajectory in a 
short sample yields the following rough rational estimate: 

QC2 = ^ = 1.111 (51) 



4.1.3 External barrier on the plasma edge 

The upper edge of the chaotic sea in Fig. (|TT|) has also been analyzed in a similar way. 
Iteration of some of the outermost points allows us to draw Fig.(|l4|) which allows us to 
select a good candidate for the most external local trajectory. The determination of a 
rational approximate for its q- value yields: 

92 

QC4 = 21 = 4,381 (52) 

For this external barrier we have successively identified on Fig.([l4|) surfaces with q — 13/3, 
35/8 and the outermost one : q — 92/21. Of course a still better precision is possible, but 
this is enough for the present purpose. 



4.1.4 Internal barrier preventing inward motion in the chaotic sea: a two-sided 
internal barrier (upper Cantorus) 

During long intermediate stages in the intermittent history of Fig. (|l3|) , the trajectory re- 
mains above what appears as an internal barrier preventing the inner motion. A similar 
analysis allows us to determine a rough rational approximate for the g-value of the inner- 
most part of this trajectory: 

qcs = ^ = 1-133 (53) 
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Figure 14: Structure of rational surfaces (islands) and chaotic zones below the edge KAM 
torus -/V(4, 2) (boundary circle) and the three convergents from below to N(i, 2): q — 13/3, 
q = 35/8, q = 92/21. This detailed drawing for L = 4.875/2vr - 0.776 represents short 
iteration times (N = 8192), performed on the most external points selected in the long 
trajectory of Figs. (11 and 13). This graph reveals different structures near the plasma 
edge, allowing to choose the most external part of the trajectory, and a rational estimate of 
the q- value of the external barrier : a KAM surface (or a robust Cantorus). 
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Figure 15: One of the last snapshot in the java animation. Bold points in this phase portrait 
(ip, 0) represent the most recent iterated points along a single magnetic line trajectory. Here, 
after filling the whole chaotic sea, the line is partly confined in the chaotic shell and then is 
sticking around the internal barrier. 



which appear to be different from the neighboring surface qci preventing upward motion 
from the inner shell. We are thus in presence of a two-sided transport barrier. It is quite 
remarkable that both sides could have been distinguished. Actually this difference clearly 
appears to the eyes on the movie jy] which has been realized on this simulation, in which we 
present a succession of snapshots to which new groups of points are added in a discontinuous 
way on each new picture : it can clearly be observed on some pictures how new, "fresh" or 
"recent" points are really aligned along a barrier which is different from the lower one. One 
of these snapshots is presented in Fig. (|l5|) . 



4.2 Expectations from theory of nonlinear dynamical systems 

From classical theory of chaos in simple nonlinear dynamical systems [|l4|, it is expected 
that the most resistant KAM torus in a q— interval of values is either the Golden number or 
at least a "noble" number (after Percival |p^7| ) . This is known to be the case in the standard 
map. When represented in a continuous fraction expansion 

[ai,a 2 ,a 3 , ....,%] = 1/ (a 2 + 1/ (a 3 + 1/... + l/a 3 )) (54) 

the Golden number has a simple [(1, )°°] coding: 

G = 1 + — ^— j- = + 1 = [1, 1, 1, 1, 1] = [(1, )°°] = 1.61803399.... (55) 

1+q 2 

and other noble numbers have a (1, )°° tail, with other integers before. Most noble numbers 
are the "most irrational", defined as those with the smallest number of integers before the 
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(1,)°° tail. 

In presence of magnetic shear the situation can be different: for instance, to the best 
of our knowledge, nobody knows why the Golden number actually plays no special role in 
the tokamap with monotonous q -profile (l6| . In other words, the special role played by the 
Golden KAM in the standard map does not seem to be conserved in other maps in presence 
of shear. We denote the most noble numbers of interest here by 

N(i,j) = [i,j,(l,)°°]=i + ^L r (56) 

3 + G 

G Z,j > 1) which represent the "next most irrational" numbers after the Golden one 
G = N(l, 1). For i = 1 we note that these noble numbers are "good milestones" in the q 
-profile since they are rather well distant, and even of measure zero JL8| , It is simple to 
prove that the values of N(l,j) are actually inserted in the q -profile between successive 
dominant island chains Q(j) = j/(J — 1) since it is easy to demonstrate that: 



Q(l+j)>N(l,j)>Q(2+j) (57) 

Better and better rational approximants to a noble number uj (in the sense of the Diophan- 
tine approximation) are obtained by truncating this infinite series at a higher and higher 
level : these rational approximants arc known to converge towards the noble number to 
and are called the "convergents", with the property "to be the closest rational to the irra- 
tional uj, compared to rationals with the same or smaller denominator" juj. The successive 
convergents to an irrational number 

uj = [at, a2,aa, ....,a,j t ...] (58) 

are defined by (see (j|J)) the series : 

[oi]=oi , [ai,a 2 ]=ai + ± , [01,03,03] = ai + a ^j_ , 



[01,02,03,04] , [01,02,03,04,05] , etc... (59) 

When looking at the most noble q— values in any given interval between go = mo /no and 
qi = rni/ni (with \mon\ — mxno\ = 1 W), it is known p3 that the most irrational number 
is 

m + G.mi 

"1 = (60) 

no + Cr.ni 

and is expected to correspond to the most robust barrier in that interval. If we now look 
at the successive intervals {Q(m + 1), Q(m + 2)} between the main rational chains of the 
dominant series q = Q(m + 1) = {m + l)/m (as observed in Fig.([ll|), we see that these 
intervals are Farey intervals [|45| covering the real axis and the most noble q— values in each 
of these intervals are precisely given by Eq. ( |56"1) : 

m+l + G(m + 2) 1 1 r , 1CC1 , , , , 

a m = m + G{7 l + l) ) ^ 1 + = 1 + m+~J = & H = N & m) (61) 



4 As a direct consequence of this fundamental recurrence relation for continued fractions, this relation 
simply expresses the fact that the two limits mo/no and mi/niof the considered interval are actually two 
successive convergents (or approximants) of some continued fraction. 



36 



Farey tree 




Rational, convergent to Noble, rational , convergent to Noble,... rational 



Figure 16: Farey tree indicating how the most noble number in a given interval between 
rationals (10/9 and 9/8, or 9/8 and 8/7) can be built with results 7V(1,8) or N(l,7). 

(where we used G + 1 = G 2 ). These most noble numbers are not frequent and simply alter- 
nate with the main rational chains. In the barrier found "around" q = 9/8 the candidates 
for robust circles are thus q = N(l,7) and 2V(1,8) since we have 8/7 > N(l,7) > 9/8 > 
N(l,8) > 10/9. 

It is worth mentioning that the convergents towards these two noble numbers are, for 
the lower Cantorus (see the lowest approximation 10/9 found in (|5l|)) : 

10 10 + 9 19 29 48 77 125 r „ n ,„ Ar/1 oN „ ^ 9 

Y-YT8--7-Y^i3^69-irr-^^^ ] = ^(1,8) = 1.116 < - 

(62) 

and for the upper Cantorus (see the lowest approximation 17/15 found in (|53|)) : 
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[1,7,(1,) 00 ]=7V(1,7) = 1.131< - (63) 



where successive approximants are alternatively from below and from above, due to the 
(1, )°° series in the Farey coding [[l5], see Fig.jT^). 

It is well known that the numerators and denominators of the successive approximants 
actually grow as a Fibonacci series = JV+i +P^, with an asymptotic growth rate 

limA^oo JV+i/fV = ^ gi ven by the golden number G. 



In Fig.(|17|), we have plotted a series of 2046 iterations along the island chains with 
rational g-values corresponding to the following approximants to the upper Cantorus (^3|) : 
q = 17/15, 43/38 and 112/99 from above, and q = 9/8, 26/23 and 69/61 from below. This 
allows us to localize the upper Cantorus on this graph as being located in the very thin 
interval located between bold dots. 
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Figure 17: The upper Cantorus for L = 4.875/27T can be localized in the central region, 
between bold dots. Here we have plotted series of 2046 iterations along each of the island 
chains with rational q values corresponding to the following approximants to the upper 
Cantorus N{1, 7) : q = 17/15, 43/38 and 112/99 from above, and q = 9/8, 26/23 and 69/61 
from below. 
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4.3 Identification of "convergents" towards these most noble Can- 
tori 

It is simple but numerically delicate to check that the barrier is composed of two Cantori 
localized precisely on these two most noble numbers : a lower Cantorus on iV(l,8) and an 
upper Cantorus on N(l, 7). Why are the most resistant KAM tori located around q = 9/8 
for this value of L could probably not be deduced on theoretical grounds. Of course, beside 
the most robust barrier described here, other, less robust barriers also exist. A detailed 
analysis of the various thin downward peaks in Fig.([l3|) reveals that the lower boundary of 
the magnetic line motion in the stochastic sea is actually oscillating in time, and does not 
always correspond to the same upper Cantorus above the island chain q — 9/8 as defined in 
Eq. (|6^) . The lower boundary of the stochastic sea can be found on different Cantori located 
most often above the island chain q = 6/5, but also on the Cantori above q — 7/6 and 8/5. 
This seems to indicate that a global barrier could rather be composed of different sticking 
regions around island chains with main rational q— values, and that these sticking regions 
are limited successively by the Cantori between them, in some kind of cascading process. 
Global motion in these boundary regions could then be described as a succession of 
approaches, like flood tide and ebb tide, allowing the magnetic line to pass through the 
successive barriers. 

During the upward motion starting from the chaotic layer, the most important barrier 
(but not the last one) nevertheless remains the upper Cantorus q = N(l,7) above the 
islands q = 9/8 : once this one has been passed, the magnetic line rapidly invades the whole 
chaotic sea. As seen in Fig.([l3|), the inverse, downward motion from the chaotic sea however 
encounters several barrier crossing, back and forth, before crossing the upper Cantorus and 
entering in the zone between the two Cantori composing the ITB described in this paper. 



4.3.1 Convergent island chains towards the two Cantori 

Localization of island chains with given q— value can be obtained numerically with great 
precision by searching hyperbolic and/or elliptic periodic points, by a numerical algorithm 
derived from a generalization of the Fletcher-Reeves method, involving the Jacobi matrix of 
the tokamap, explained in Appendix B. Localizing the position of a noble Cantorus can 
only be achieved as a limiting procedure, by localizing the series of its convergents, as given 



by Eqs. fcg, 63) 



In order to check the above predictions we have sampled (in a very long trajectory of 2.10 9 
iterations) various sections of trajectory (a) reaching the highest tp— values in the chaotic 
shell (below the lower Cantorus), and (b) reaching the lowest ip— values in the chaotic sea 
(above the upper Cantorus): such trajectories remain indeed well separated, by the width 
of the barrier. Then we have checked that any of these sections of trajectory (a) remains 
indeed below the limit of convergence of iV(l, 8) of the lower Cantorus, more precisely below 
the limit of the convergents from below : 



10 29 77 202 529 1385 
9 ' 26 ' 69 ' 181 ' 474 ' 1241 > 



iV(l,8) (from below) (64) 



and that any of these sections of trajectory (b) remains indeed above the limit of convergence 
of N(l,7) of the upper Cantorus, more precisely above the limit of the convergents from 
above : 

17 43 112 293 767 2008 =± N(l 7) (f rom above) (65) 

15' 38' 99 ' 259' 678' 1775 — ^ iV V 1 ' ') l uuul auovcj yuo j 

In other words we have checked that observed trajectories in the chaotic shell and in the 
chaotic sea remain actually always on their own side of the pair of Cantori, which localizes 
the barrier. 
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4.3.2 Convergent island chains towards noble KAM's on the edge and around 
the plasma core 

By the same method we have also identified the q— value of the robust boundary circle 
forming the plasma edge as being equal to TV (4, 2) which appears to be the most irrational 
between q = 4 and q = 4.5, according to Eq. (pG) : 

4 4 + 9 _ 13 22 35 57 92 149 

T 1 + 2 ~T^T^T^13^21^ ~34~' "' ^ 

[4,2,(l,) 0o ]=iV(4,2) = 4.381<^ (66) 

in agreement with (|52|). There should probably exist a stronger barrier out of the "plasma 
edge" on iV(4, 1) ~ 4.618, but this one remains out of reach for trajectories starting from 
the plasma bulk at the considered value L — 4.875/27T since the KAM on the edge at 
q = N(4.2) < N(4, 1) has not yet been broken. 

We have also identified the q— value of the robust boundary circle protecting the plasma 
core as being equal to N(l, 11), the most irrational between q = 1 and q = 12/11, according 
toEq. @: 

1 1 + 12 13 25 38 63 101 



1 1 + 11 12 23 35 58 93 ' 

=> [1, 11, (1, )°°] = N(l, 11) = 1.086 < H (67) 

which is in poor agreement with the above rapid estimate d50j), but correct to the first three 
digits. 

4.3.3 The ITB: a double sided barrier around q = 9/8 

The final scheme which results from the above measurements on a very long tokamap tra- 
jectory at L — 4.875/27T is the following. A magnetic line with an initial condition inside the 
inner shell (IS) actually has an inward motion limited by a robust KAM torus protecting 
the plasma core (or a Cantorus ?) at 

q IS > N(l, 11) = 1.086 > H (68) 

and an outward motion limited by a semi-permeable Cantorus at: 

q IS < N(l, 8) = 1.116 (69) 

This numerical analysis shows that noble Cantori are good candidates to be identified with 
internal transport barriers, as could have been anticipated from the relation between KAM 
theory and number theory [jisj. 

Once arrived in the main chaotic sea (CS) extending up to the plasma edge, the magnetic 
line wanders around island remnants, remains stuck around some of them (here in Figs.([l3|, 
fjr| ) around the q = 5/2 chain), but its inward motion is limited by a semi-permeable 
Cantorus 

q CS > N{1, 7) = 1.131 > | (70) 
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The measurements indicate that this inward-motion limit in the chaotic sea is different from 
(and located above) the outward motion limit of the IS : 

N(l, 8) = 1.116 < I = 1.125 < iV(l, 7) = 1.131 (71) 
8 

We have thus observed the existence of a two-sided transport barrier around q — 9/8 and 
limited on each side by a Cantorus, respectively at q = N(l, 8) and q = N(l, 7). 

On the other end, the outward motion in the chaotic sea is limited by a curve which will 
be proved in Section (||) to be a KAM surface at 

qcs < N{4, 2) = 4.382 (72) 

which is the most irrational number between 4 and 4.5. This external KAM surface thus 
appears indeed a good candidate for the external barrier observed here. 



We are thus in presence of a double-sided transport barrier, composed of two Cantori, the 
lower Cantorus with noble q— value JV(1,8) = 1.116, the upper Cantorus with noble value 
N(l,7) = 1.131. These Cantori appear on Fig.©. From @ we have N(l,7) > Q(9) > 
iV(l, 8) , which shows that the rational surface q — Q(9) — 9/8 is actually between these two 
Cantori and thus inside the transport barrier. In this sense, one can say that this internal 
transport barrier is actually located around a dominant rational surface, in spite of the fact 
that it is actually composed of two irrational, noble Cantori. 



4.3.4 Experimental localization of barriers in the Tokamak q -profile 

From the experimental point of view, transport barriers in tokamaks are indeed generally 
observed "around" the main rational q- values [|). 

In the RTP tokamak, with a wide range of g-values (0.8 =>■ 5) in a reversed sheared 
profile, it has been observed, by varying the heat deposition radius p dep of off-axis ECH 
heating, that the central electron temperature 2^(0) decreases by a series of plateaux H. It 
is observed that the values of q m of the different plateaux fall in half-integer bands, i.e. q m 
crosses a half-integer value each time the discharge transits from one plateau to the next. 
Since these transitions correspond to the loss of a TB, these authors deduced that "the 
barriers are associated with half-integer values of q ." 

It has been reported that some " ears" appear in the T e profile (appearance of one bump 
on each side of the central value) when the heat deposition radius p dep is localized inside 
a transition between two plateaux M . We note that the existence of such " ears" could be 
an indication in favor of the existence of a double-sided semi-permeable TB around p dep , 



as found in the previous Section (4.3.3). Such "ears" appear to be unstable and to crash 
in a repetitive fashion, showing that the barrier can indeed be crossed and appears to be 
permeable. 

The central sawteeth allows to place the first barrier "near" q = 1 ; off-axis sawteeth 
indicate other barriers near q = 3/2, 2 and 3. Remaining barriers are attributed to q — 4/3 
and 5/2. 

In all cases the barriers observed in experiments are associated with those "dominant 
rational" g-values, but a specific experimental work could hardly have been done in order to 
determine more precisely the possible role of noble or other irrational values around these 
" dominant rationals" . 

Up to now it has been generally admitted 0], ||, JhJ that internal transport barriers 
could correspond to rational g-values, where primary island appear. Chaotic motion can be 
observed between these primary rational islands but mainly around the hyperbolic points. 
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The picture we obtain here is rather different. It is well known that, in presence of several 
island chains, secondary islands appear and accumulate around irrational surfaces where 
KAM surfaces finally broke themselves into discrete pieces forming a Cantorus. This is 
precisely the location where we have found internal barriers: they appear as Cantori on 
irrational surfaces rather than rational surfaces. 

This finding can in turn be helpful to build transport models like the g-comb model p3| , 
|24j. From experimental considerations, such models define indeed barriers are localized 
around main rational g-values, but the width if the barriers has still to be defined. We 
propose to localize the edges of each barrier on the most irrational value in the prescribed 
domain. 

4.4 Last barrier in the standard map: the golden KAM 

In order to illustrate the destruction of a robust barrier, there exists one example in which 
a KAM surface can be represented just at the critical point where it is broken into a Can- 
torus. This is the case of the standard map [^7j where the critical value of the stochasticity 
parameter is known with a very high precision. It is interesting to remind that the last 
KAM surface in that case has been identified to be the Golden KAM with q- value equal to 
G, the golden number. The breaking of this surface and its transformation into a Cantorus 
is known to occur at a critical value of the stochasticity parameter which is known to 
be given by K — Kc — 0.971635406. Because such a drawing is not easily found in the 
literature for this value of K, it appears worthwhile to present the trajectory following this 
critical curve. In Fig.([l8|) we present part of the phase portrait of the standard map for 
K = Kc- Two series of islands remnants q = 8/5 and q — 5/3 can be seen, corresponding 
to two convergents towards the Golden number G = lim of (|, |, | , ...) . The last KAM 
with q = G is exhibited inbetween, and it appears as the last existing non-chaotized KAM 
surface, surrounded by chaotic layers around rational islands. 

This example shows the detailed structure of an irrational KAM barrier surrounded by 
chaotic, permeable zones. 



4.5 Perturbed q -profile in the tokamap 
4.5.1 Exact perturbed ^-profile 

In order to understand the radial positions of these barriers on the exact q -profile of the 
tokamap, we have succeeded to draw a one-dimensional g-profile of the magnetic surfaces, not 
as function of the radius (since perturbed surfaces are not circular anymore), but as function 
of the distance X between the polar axis and their intersection point in the equatorial plane. 
We have calculated the intersections of several trajectories with the equatorial plane (at 9 = 
and 9 — 0.5) and measured their g-value, by computing the average increase of the poloidal 
angle according to Eq.j4§|). 

In this way we obtain a rather precise profile of the perturbed q- values of various magnetic 
surfaces or island chains, represented at the two (non-symmetrical) points where they cross 
the equatorial plane, see Fig.(|l9|). Let us recall that the radius r represented on the various 
figures is defined by r = (see Eq.(|36|). This variable r varies from to a/2 on the edge 
of the plasma, while the variable x — \p4> = r/^/2 varies from to 1 and represents the 
reduced radial coordinate with respect to the small radius of the torus. 




(73) 
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Figure 18: Part of the phase portrait of the standard map for K — Kc- Two series of islands 
q = 8/5 and q = 5/3 can be seen, corresponding to two convergents towards the Golden 
number G. The last KAM with q = G is exhibited inbetween, and it appears as the last 
existing non-chaotized KAM surface; its initial conditions are / = 0.52969, 8 — 0.75. 
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Figure 19: Calculated points (diamonds) of the perturbed g-profile (around the magnetic 
axis) in the equatorial plane for L = 4.875/2-7T ~ 0.776, with the plasma extending from 
X = — 1 to +1, showing overall agreement with the unperturbed profile (dotted line), except 
for the the appearance of a local maximum on the magnetic axis (XMA). 



The abscissa in Fig.([l9|) is represented between X = x cosO = 1 on the weak field side and 
X = — 1 on the strong field side, as compared with the unperturbed imposed q -profile (Eq. 



The measured points roughly agree with the unperturbed curve (Fig. (19)). But a detailed 
drawing in Fig.(^) reveals a rather unexpected point: the perturbed safety factor profile 
in the equatorial plane exhibits a local maximum on the magnetic axis, and a minimum on 
both sides. 

The measured g-value near the magnetic axis appears to be 

qma =s 1.1346 ^ (74) 
15 

This non-monotonous perturbed q -profile , which is spontaneously created by the magnetic 
perturbation, could be the reason for the appearance of transport barriers in the tokamap. 



4.5.2 An ITB in the rev-tokamap 

Transport barriers are very frequent in reversed shear situations, and they also appear in 
the tokamap model for such cases. We recall indeed that a transport barrier has already 
been shown to occur around the minimum of a reversed shear q -profile (" revtokamap" p7|). 
In the discussion of Fig. 4 in Ref. p7[ , for L = 2.8/27T, it is observed that "the chaotic 
region is sharply bounded from below by a KAM barrier" located at 1/q = 0.5772. We 
note that this value is only 0.5% away from a noble value 1/q = 0.5802 which corresponds 
to q = [1,1,2,(1,)°°] = 1 + (1/(1 + 1/(2 + 1/G))) = 1.7236068... which is larger than the 
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Figure 20: Detail of the perturbed g-profile (measured around the magnetic axis) in the 
equatorial plane (circles) for L = 4.875/27T ~ 0.776, as compared with the unperturbed 
profile (dotted line) , showing the appearance of a local maximum on the magnetic axis, and 
of two minima. Horizontal levels represent (from top to botom) the g-values N(l, 7) = 1.131 
of the upper Cantorus (which corresponds approximately to the higher q- values of the bump) , 
N(l, 8) = 1.116 of the lower Cantorus, and N(l, 11) = 1.086 of the KAM barrier protecting 
the plasma core, along with the extension of these curves with the equatorial axis X. The 
value on the magnetic axis is Qma = 1-1346 > iV(l,7) and lies at the center of a profile 
which is locally hollow at very small distances from the magnetic axis (not visible at this 
scale) 
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Golden number and appears from (^(J) to be the most irrational number between 5/3 and 
7/4. Within a precision of 0.5%, the transport barrier found in |37| could thus well be 
located on a noble value again, even in the previous situation of a reversed magnetic shear 
profile. 



4.5.3 Inverse shear in the tokamap 

The inner shell, which has been observed between qci — 1.080 and qc2 = 1.111, actually 
involves q- values lower than that measured on the magnetic axis (Eq. (|74|)), down to q ~ 1.07 
(see Fig.(p0[), precisely because of the non-monotonous profile. 

One can see that the g-values decrease from the magnetic axis towards a separatrix 
passing through the polar axis (on which two hyperbolic points are superposed, as long as 
L < 1). This whole region inside the separatrix would appear however to have q = 1 if 
and only if q was computed around the polar (geometrical) axis, but this way of counting 
has not much physical meaning, besides showing that this plasma core could appear as an 
island m = n — 1 in the reference frame of the polar axis. As a result, trajectories inside 
the chaotic shell have actually q— values (around the magnetic axis) which decrease from the 
boundary circle for those trajectories which do not enclose the polar axis (ip = 0,8 = 0), 
but which increase for those which encircle the polar axis, up to the value N(l, 8) of the 
lower Cantorus. In the chaotic sea, the q -profile is monotonously growing : trajectories 
have q— values higher than the upper Cantorus and lower than the boundary circle on the 
edge, which is the expected situation. 

In summary, starting from the magnetic axis, we first find a regular zone, with a decreas- 
ing perturbed q -profile away from the magnetic axis (where quA — 1.13461), then a robust 
boundary circle located on q = N(l, 11) protecting the regular plasma core, then decreasing 
q— values up to the separatrix, then a regular increase of the q -profile with a q = 10/9 
rational chain, a lower semi-permeable lower Cantorus on q = iV(l,8), the island remnants 
q = 9/8, then the upper Cantorus on q = N(l, 7) below the rational chain q = 8/7, etc... 
up to the robust KAM torus on the edge at N(A, 2). The two robust barriers and the two 
Cantori are represented on Fig.(|2l|) along with a part of each chaotic region. 



4.6 Sticking measurements 

On the other hand, we have seen in Fig. (|l3|) that a long trapping stage is observed in the 
chaotic zone surrounding the island remnants q — 5/2. For a much longer trajectory (10 11 
iterations) one observes that such temporary but long trapping can occur around almost all 
rational and specially around the main islands known to play a role in the tokamaks: q = 2/1, 
5/2 and 3/1. In Fig.(p2]) we have computed the normalized histogram showing the occurrence 
H of sticking times longer than 10 5 (in a run of 10 11 iterations with L = 4.875/27r), obtained 
by a sliding average, as a function of the radial position if) of the visited island chain. The 
values of ip are taken within small intervals of 2.10 -4 . This graph exhibits the frequent 
occurrence of long sticking times in and around the main rational chains : the two bumps 
correspond to widespread ip values of a line wandering in the chaotic shell (around ip ~ 0.1) 
and in the chaotic sea (around ip ~ 0.4 to 0.8). These are the main two chaotic zones. On the 
other hand one remarks sharp peaks representing long sticking events at the corresponding 
ip values. In the chaotic sea bump, we identified very precisely (by measuring the exact 
average-^ value of the corresponding q) long and frequent sticking events around island 
remnants q — 2, 5/2 and 3 (with satellite peaks at neighboring rationals), etc... This kind 
of graph yields a precise measurement of the richness of trapping phenomena. 
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radius 




angle 

Figure 21: The main noble transport barriers are represented in bold lines in normalized co- 
ordinates: "radius"= x = p/a = y 7 ^ and "angle"= 9. The robust KAM N(l, 11) separating 
the central protected plasma core (in white) from the chaotic shell. The two semi-permeable 
Cantori JV(1, 8) and N(l, 7) form an internal barrier resulting in a very slow and intermit- 
tent motion towards the chaotic sea. The robust KAM torus on the plasma edge N(4, 2) has 
been identified to have a vanishing flux. 
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Figure 22: Normalized histogram showing the occurrence H of sticking times longer than 10 5 
(in a run of 10 11 iterations with L = 4.875/27r), obtained by a sliding average, as function 
of the radial position ip of the visited island chain. The values of ip are taken within small 
intervals of 2.10~ 4 . This graph exhibits the frequent occurrence of long sticking times in 
and around the main rational chains : the two bumps correspond to widespread tp values 
of a line wandering in the main two chaotic zones: the chaotic shell (around ip ~ 0.1) and 
in the chaotic sea (around ip ~ 0.4 to 0.8). 
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It would be interesting to study the characteristics of the other secondary barriers, sur- 
rounding the " holes" in the phase portrait Fig. ( |i"l| ) and to determine if such local barriers 
also correspond to noble Cantori in the profile of the local q- value (around the local magnetic 
axis in the center of the island remnant). For instance a sticking region is easy to observe 
along the edge on the q — 9/8 islands. 

The interest of this new description of ITB's in tokamaks is that mathematical tools 
exist for the description of the motion across such Cantori ("turnstiles" etc... 46 ), and 
these tools can be used to study transport across internal transport barriers in confined 
plasmas. Calculation of the flux through a Cantorus, defined by the limit of its convergents, 
can be performed by using the turnstile mechanism ^6|. This subject is presented in the 
next Section ^|. 



5 Calculation of the flux across noble barriers 

The study of a single long orbit points out some transport phenomena which occur between 
some phase-space's zones separated by internal transport barriers (ITB), but it does not 
give information about the flux through these barriers ( i.e. the area of the set formed by 
points which pass through a ITB at every iteration). This information will be obtained 
in this section from Mather's theorem. Using Greene's conjecture we confirm the observed 
existence, for L = 4, ^J 5 ~ 0.776, of the invariant circle on the plasma edge having the 
rotation number u>i = N ^ 2 ) an< ^ °^ a Cantorus having the rotation number u>2 — jy(i,7) • 
Both Greene's and Mather's approaches use the intimate connection between quasiperiodic 
and periodic orbits. 

5.1 Periodic and quasiperiodic orbits. The action principle 

Let us remind some definitions and useful basic results in the theory of dynamical systems. 

The" lift" of the map T : [0, 1) X R ->[0, 1) x R, T (9, ip) = (9' (mod 1) , ip') is the appli- 
cation T:RxR^R defined by T (9, ip) = ', if}') (because the (mod 1) is not applied for 
9' , the lift T allows to take into account the number of turns along 9). 

The map T : [0, 1) x R ->[0, 1) x R is called a twist map if ^ for all (9, t/>) G 

[0, 1) x R. It is a left (respectively right) twist map if " < (respectively 96 > 0) 
for all (9, ip) G [0, 1) x R, which means that the perturbed winding number is a monotonously 
decreasing (respectively increasing) function of ip. 

A point (6*o, ipo) is a periodic point of type (n, m), with periodicity to, if (9o,ip ) 
= (9q + n,ip ). The rotation number of a periodic orbit exists and is equal to l/2tt = n/m 



(i.e. the limit in Eq.(48) exists); it is clearly independent of the initial point on a periodic 
chain. This global property persists for invariant curve, under some restriction on the map. 
The rotation number may also be irrational in this case. Beside periodic orbits and invariant 
curves (called invariant circles, for topological reasons), an intermediate case of invariant 
set exists, on which rotation number exists and is irrational: the Cantori. The dynamics 
on these three types of invariant set is generally quasiperiodic: the "time" dependence can 
be expressed by a generalized Fourier series, with rational frequencies for rational u, and 
irrational incommensurate frequencies in the remaining cases (see Ref. |l7j or Rcf. |Q for 
details). Both periodic and quasiperiodic points can be obtained as stationary points in the 
action principle. 

The " action principle" for maps is the analogous of the Lagrangian variational principle 
in continuous dynamics. In the case of Hamiltonian twist maps, the action generating 



function F a (defined in Section (2.2)) plays the role of a Lagrangian for discrete systems. 
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Theorem 1 (action principle for periodic orbits) Let T : S 1 x R — > S 1 x R be an area pre- 
serving twist map , F a : RxR — > R, its action generating function and {(8q, ipy) , {8^ n _ 1 , Vv, 
a periodic orbit of T of period m. Then {8q, ...,8^ n _ 1 ^ is a stationary point of the action 

A (8 , 6> m _i) — F a (6 , 81) + F a (8 1 ,8 2 ) + ... + F a (6> m _i, O ) (75) 

A very simple proof is presented in Ref. (4^] 

In 1927 Birkhoff |30) showed that every area-preserving twist map has at least two 
periodical orbits of type (n, m) for each rational winding number - = j n an appropriate 
interval (called the twist interval). For area-preserving twist map it can be proved (see Ref. 
|47[| p. 38 for commentaries) that, for each rational ^ in the twist interval described by the 
map, there exists at least one periodic orbit of type (n, m) which extremizes A (it is called 
extremizing orbit) and at least one periodic orbit of type (n, to) which is a saddle point of 
A (it is called a maxmin orbit). 

In order to study the linear stability properties of a (n, to)- periodic orbit one computes 
the multipliers of the orbit {i.e. the eigenvalues A and j of the Jacobi matrix associated to 
T o T o .... o T in an arbitrary point of the orbit) or the residue 

m — times 

2--A_-_l/A (76) 

If R < (i.e. A > , j > 0) the orbit is formed by direct hyperbolic points (which are 
unstable). If R 6 (0, 1) (i.e. A and j are complex conjugate numbers and |A| = |^| =1) the 
orbit is formed by elliptic points (which are stable). If R > 1 (i.e. A < 0, j < 0) the orbit 
is formed by inverse hyperbolic points (which are unstable). 



In Ref. |3l|] it was proved that every extremizing orbit has negative residue and that 
every maxmin orbit has a positive residue. It results that every extremizing orbit is formed 
by direct hyperbolic points and the maxmin orbit is formed by elliptic points (if R £ (0, 1)) 
or by inverse hyperbolic points (if R > 1). 



The action principle for quasiperiodic orbits was proposed by Percival (in Ref. 52 | V In 
1920 Birkhoff proved that any quasiperiodic orbit is the graph of a function (see Ref. f53|). 
It can be written in the form 

{(«(*),#)), teR} (77) 

where 8 : R — * R is an increasing function having the periodicity property 6 (t + 1) = 6 (t)+l 
for all tell, and where 

f)F 

m = -^f[e(t),0(t + ")} (78) 

In these terms the action principle can be written as follows. 

Theorem 2 ( the action principle for quasiperiodic orbits) Let T : S 1 x R — > S 1 x R be 
an area preserving twist map, F a : R 1 x R 1 — > R its action generating function, to € (0, 1) 
an irrational number and {(8 (t) , ip (t)) / t G R} a quasiperiodic orbit having the rotation 
number u). Then 8 is a stationary point of the functional 

1 

A (8) = [F a [8 (t) , 8(t) + lu] dt (79) 
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For twist maps, Mather (see Refs. and p5| ) proved the existence of a stationary and 
moreover " extremizing" function 9 ext which extremizes the functional A (|79|). It gives rise 
to an invariant set (C). Its equation can be obtained from ([77|, |7^). If 6* e xt is continuous then 
(C) is an invariant circle. If 9 ey ± is not continuous (but has a countable set of discontinuities 
because it is monotonous), then the closure (the set of the limit points) of (C) would be a 
Cantor set, called "Cantorus". A very good survey on this problem is given in Ref. Q. 

5.2 Study of quasiperiodic orbits via periodic orbits 

The invariant set (invariant circle or Cantorus) having the irrational rotation number lu 
was obtained using a stationary point 9 ext of the action functional A (6) defined in ([79]), 
but it is also the limit circle of a sequence of periodic orbits having rotation numbers — 
, when lim — = lu. High order periodic orbits may thus be considered as good enough 

approximations for the invariant set and the study of their properties gives information 
about the invariant set properties. So, an irrational magnetic surface can be described as 
the limit of its rational convergents by observing higher and higher order periodic motions. 

There are at least two approaches to study the connection between these invariant sets 
and the periodic orbits: the Greene's conjecture (which relates the existence of the invariant 
set to the stability of a particular sequence of periodic orbits) and the Mather's theorem 
(which gives necessary and sufficient conditions for the existence of the invariant set and 
computes the flux through it). 

For a noble number lu we will denote by {j^j the sequence of its convergents (see 



Sections 4.2 and 4.3 for definitions and computations). We denote by i?+ the residue of a 
maxmin orbit (passing through elliptic or inverse hyperbolic points) and R~ the residue of 
an extremizing orbit (passing through direct hyperbolic points) of type {n v ,m v ). 

Greene's conjecture (in Refs. [|5(| and J57J) predicts that there would be three situations: 

- subcritical ( — > and R~ — > ); in this case there is sequence of island chains 
converging to a smooth invariant circle having the rotation number lu. 

- critical ( i?+ and R~ -/> but they are bounded in (0, oo) and in (— oo,0)), 
respectively; in this case there is a sequence of island chains converging to a non smooth 
invariant circle having the rotation number lu. 

- supercritical ( i?+ — > oo and R~ — ► — oo); in this case there is no invariant circle having 
the rotation number lu, but there is a Cantorus with this rotation number lu. 

In Ref. |5a] a partial proof is presented, but the conjecture is not yet completely proved. 



On the other hand, in Ref. (59[ Mather gives an equivalent condition to the existence of 
an invariant circle having the rotation number lu. For a left twist map and a rational ^ in 
the twist interval one can define 

^•^-n/rn ^extremizing -A maxmin 

(80) 

the difference of action A between the extremizing orbit and the maxmin orbit of type 
(n, m) . In ( |S0"| ) | AAn/m | is the escape area from a domain bounded by an arbitrary curve 
(C) passing by the m points of the maxmin orbit and by the m points of the extremizing 
orbit of type (n, m) .Let us remind that the escape area (the flux) from the domain Rm) 
bounded by the closed curve (C), under the map T, is defined by 

L(R Cl T) d ^ Area(T(R c )-Rc) (81) 
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T <V R c 




Figure 23: The figure visualizes the escape area from a bounded domain. The curve T(C) 
is the image of the closed curve (C) after one iteration by the map T . The domain Rc is 
the domain bounded by (C). The set T(C) fl (C) consists in three unstable periodic points 
(namely U(l), U(2) and U(3)) and in three stable points (namely 5(1), 5(2) and 5(3)). 



In Ref.[|46| a very interesting result is proved for right twist maps. It can be easily 
adapted for the twist maps as follows. If (C) n T (C) consists in 2m points of an extremizing 
orbit and of a maxmin orbit of type (n, m) then 



L(R C ,T) = 

I Axtrcmizing ~ A maxmin 



m—1 

3=0 



j ' 



1 J + 



\) F a (t 



3 ,C J+1 



(82) 
(83) 



| ^A n j jyi | 

In this formula, \8 < f'\'4 > ^ > ) an( i (^•^iV'j^)) f° r 3 e {1> 2, m}, are the points of the 
extremizing orbit and (respectively) of the maxmin orbit of type (n, m). An example with 
periodicity m — 3 is illustrated in Fig (23). 



Mather's theorem states that there is an invariant circle having the rotation number u> 
if and only if lim AA n i m = 0. Otherwise, the sequence (A^4 n / m ) converges to a 
nonzero limit A_A W . In this latter case there is no invariant circle of rotation number lo, but 
there is a Cantorus with this rotation number, and |A^4 W | is the flux through the Cantorus. 
It results that the flux through a Cantorus can be approximated by | AAn v / mv \ for large 
enough values of v. 



The system ( |l6| ) gives rise to an unique area preserving map T : R x R + — > R x R + 
such that T(0„, = V'jz+i) • If 7 1 is a twist map and if the analytical form of F a can 
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be obtained, the formula ( |S2| ) can be used in order to compute AA n ^/ mv (as in the case of 
the standard map) . But there are a lot of systems (the tokamap- model for example) whose 
analytical form of F a can not be derived. 



In order to obtain L(Rc,T), a direct computation can be done. If we denote by 
/ (VV+i> ^f) = Fo (ipv+i) + L5F (ip u+ i,O v ) , the generating function flTij) becomes 

f (v>„+i. e v ) = i) v+l -e v + f ty v +i,o«) (84) 

and the system ( |l6| ) can be written as 



One observes that / (Vv+ij ^) in (85) plays the role of a discrete Hamiltonian for the map, 
which is periodic in Q v . This is not the case for F (Vv+ij ®v) in (gj) - where i* 1 is the explicit 
generating function of the tokamap, given in (|30|). 



By repeating the arguments applied to the calculation of the escape area in terms of 
the action function, we deduce the following expression in terms of the periodic function 
/ (Vv+i A): 

= e • - / (86) 

j=o \ uip v+1 ) 



j'=o 



Here f^"^)^"^ an d ffl^'^M, for j £ {1,2,..., to}, are the points of the extremizing 
orbit and (respectively) of the maxmin orbit of type (n, to) contained in C . 



Because of the relation: 



: Vi (^+i - *i) - F (^i+i. e i) + Vi • B i z 
= • 6 j+1 - F 6j) = F a {piJj+i) 



it can easily be shown that the formula ( p6[ ) involving only periodic functions is actually 
equivalent to the following form: 

m— 1 

where F is the (non periodic) explicit generating function of the tokamap, given in (|3^). In 
this formula (p8|), ^^iV^"') an< 4 (^j^' ' ^ or e {1,2,..., ?n}, are the points of the 
extremizing orbit and (respectively) of the maxmin orbit of type (n,m). 



53 



The new formula (|36|) for the escape area can be used when only the mixed generating 
function is given. It works also in the case of the non-twist maps, when the action generated 
function involved in ( |8^ ) cannot be defined and the formula (82) can not be used. The 
representation of mixed generating function F given by (|l4|) is recommended because / is 
periodic in 8, unlike F and F a . This explicit periodicity prevents programing mistakes and 
increases the numerical stability. 



5.3 Results about the tokamap, a twist map 

The discrete system (pl|)-(p4]) is generated by the original "tokamap" namely T : [0, 1) 
R+ -> [0,1) x R+, defined by T(6,ip) = (0',ip') 



i/i-l-Lsm (2ir0) + J^j} - 1 - L sin {2ir0)f + Axp 



9' = 8 + \ (2 - V') 2 - 2ip' + (V') 



L cos(2tt0) 

2¥ ' JT+^V 



(mod 1) 



(89) 



It can be checked that the tokamap is a left twist map. Indeed, a simple derivation yields 



98' (fl.VQ 
dip 



dip 



5» <0 



(90) 



Numerical computations show that the first bracket in ( p0| ) is negative in [0, 1) x [0, 1] 
which is the interesting zone in the phase space. On the other hand 



dip' 
7hp 



ip + l-L sin (2tt0) 



\J(ip - 1 - L sin (2tt6»)) 2 + Aip 
for all (8, ip) £ [0, 1) x [0, 1] and L < 1. As a result, the inequality 



1 

> 2 



(91) 



is obtained. 



The analytical form of the action generating function F a (8 u ,8 u+ i) defined in (|13| ) for 
the Tokamap cannot be derived analytically from the map equations (^9|). To obtain it 
one has to solve the equation ( |32"| ) for the unknown Vv+i; using the g-profile ( [l5| , |39| ), and 
to substitute it in (|l3|). It can be proved that the equation ( |32| ) has an unique solution 
ip u+ i = ip v -Li {8u,9v+i) but this solution cannot be obtained through analytical methods. 
However, numerical values of F a (9 V , 9 v +\) can be computed. In order to compute AA n / m , 
the formula (|8^) will be used. 



5.4 Numerical results for computing the fluxes in the case L = 

4.875/2tt 

The numerical results presented in this paragraph actually confirm the existence of a KAM 
barrier with rotation number u>i = l/iV(4, 2), as noticed above. We also evaluate the flux 
through the upper Cantorus of the internal transport barrier, which has the rotation number 

L0 2 = 1/A(1,7). 

In computing these fluxes the identification of long periodic orbits is important. Since 
the phase portrait of the tokamap has no symmetries it is quite difficult to localize long 
periodical orbits. This can however be done by using an algorithm, described in Appendix 
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B, which is based on the Fletcher-Reeves method. In spite of the large number of digits 
used in the calculations (18 digits), the roundoff errors (see a discussion p. 276 in |19[|) still 
introduce numerical deviations and errors in computing long periodic orbits, for the well- 
known reason that hyperbolic periodic points are unstable and the system exhibits a strongly 
sensitivity to the initial conditions. 

For the (observed) KAM barrier having the rotation number u>i = l/N (4, 2) the sequence 
of the convergents was obtained in (^). The results of the numerical computations for i?+ , 
\R~\ and AA n ^/ mv are presented in Table (|92j). 



// 




R+ 


R v 




^•^-n^ jra v 




1 


5/22 


0.2321833 


-0 


2369255 


1.7174326e- 


05 


2 


8/35 


0.2571793 


-0 


2626238 


4.4808639e- 


00 


3 


13/57 


0.2376864 


-0 


2425847 


9.5294768e- 


07 


4 


21/92 


0.2446923 


-0 


2497698 


2.2566047e- 


07 


5 


34/149 


0.2323534 


-0 


2370065 


4.9617970e- 


08 


G 


55/241 


0.2271472 


-0 


2315720 


1.1173445e- 


08 


7 


89/390 


0.2146109 


-0 


2146109 


2.4016014e- 


09 


8 


144/631 


0.1910067 


-0 


0612393 






9 


233/1021 


0.1604416 


-0 


1627059 






10 


377/1652 


0.1219058 


-0 


1232303 






11 


610/2673 


0.0776218 


-0 


0781681 







The periodic points were computed with an error of the order 10 13 . We can notice that 
R+ is decreasing from R+ = 0.2321833 to i?+ = 0.077621821, and that R~ is increasing 
from R^ = -0.326925525 to i? u = -0.078168125 in 11 steps, approaching to 0. The 
Greene's conjecture (subcritical case) is fulfilled. Probably because L = 4.875/2-7T is near 



the threshold value for breaking up the KAM barrier (see Section 3.1.6 and Conjecture in 
Section [T^) the convergence is very slow. 

We observe that the sequence of values for the fluxes (AA nv / mi/ ) v& ^s is rapidly decreas- 
ing, hopefully towards 0, and this allows us to argue that the flux through this noble KAM 
surface q = N(4, 2) is indeed zero, indicating the existence of a strongly resistant barrier 
on the edge of the plasma for L = 4.875/2-7T, as observed in the simulations. 

The computer errors enable us to evaluate correctly the fluxes through a curve passing 
through the points of a Poincare-Birkhoff chain containing more than 2x390 = 780 points, 
even the formula for the escape area can be used. But the error in computing the fluxes 
increases from 1.862 -10" 15 (for AA niiy/mi ) to 2.081 TO -13 (for AA ni / m7 ) and becomes of the 
same order of magnitude as the flux for v > 7. From this high periodicity, the computations 
are not precise. 

For the upper Cantorus having the rotation number ui^ — j^k 7 ) , the sequence of con- 



vergents was obtained in ([33]). The results of the numerical computations for \R V \ and 
A-4„^/ mv are presented in Table (p3[). 
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V 


/ 

n v jm v 


r>4- 

R+ 


R v 




* 10 5 


1 

1 


s/y 


U.Uz I zUoy40D 


fl 07C01 7/1 1 a 

— U.zVozl<410 


4. ( obUe — 


Uo 


2 


15/17 


0.3082029734 


- 0.316420328 


7.9925e - 


06 


3 


23/26 


0.3646525501 


- 0.374840518 


2.4865e - 


06 


4 


38/43 


0.4459480918 


- 0.462159188 


6.5746e- 


07 


5 


61/69 


0.6621904100 


- 0.694999295 


2.2083e - 


07 


G 


99/112 


1.2004032200 


- 1.297669846 


8.6218e- 


08 


7 


160/181 


3.3011913350 


- 3.814698002 


4.6269e - 


08 


8 


259/293 


17.1200315801 


- 22.284309839 


3.4357e - 


08 


9 


419/474 


221.306197153 


- 166.19991214 






10 


678/767 


1066.799716960 


- 546.73330000 






11 


1097/1241 


1971.040117100 


-1368.80947500 
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We notice that Greene's conjecture (the supercritical case) again is fulfilled: increases 
quickly (hopefully to +oo) and R~ decreases rapidly (hopefully to -co). Due to the com- 
puter's errors, we can not evaluate the fluxes through a curve passing through the points of 
a Poincare-Birkhoff chain containing more than 2x293 = 586 points. 

Because the escape areas' sequence is decreasing with a decreasing slope, we may 
consider that it goes to a constant, and that its limit is less than 3.4357 • 10 -8 and is of the 
order of 10~ 8 . This is the value of the magnetic field lines' flux through the internal barrier 
having the rotation number iV(l,7). It means that, at every iteration, the magnetic field 
lines passing by points which do occupy a surface (in the Poincare section) having an area 
of approximately 3.4357 • 10~ 8 pass through the Cantorus, in the chaotic sea. In the same 
time some other magnetic lines (located in a set with the same area in the chaotic sea) come 
inside the region bounded by the Cantorus. 

We can also observe that we can find elliptic points in every neighborhood of the KAM 
barrier, but near the partial transport barrier the high order convergents are direct or 
inverse hyperbolic, with exponentially increasing (with the periodicity order) residues. It 
results from this that the partial transport barrier is located on a more chaotic zone, as 
compared to the KAM barrier, where there are many zones with regular dynamics (near the 
elliptic points) which are embedded in the chaotic zone. The chaoticity of the chaotic zone 
is reduced in the vicinity of the KAM, because the residues of the direct hyperbolic points 
approach to zero. 

The results can be summarized as follows. For the KAM barrier, the positive and the 
negative residues converges slowly to 0, so the subcritical case in the Greene's conjecture 
is verified . The possible explanation for this slow convergence is that the value of the 
stochasticity parameter K = 4.875 is very near the threshold for large scale dispersion due 
to the breaking up the edge barrier. By computing the periodic orbit with an error of 10 3 , 
we obtain still decreasing values of the flux of the order of 10 -9 , actually smaller than the 
flux across the Cantorus. So, we may expect the convergence towards 0. 

For the Cantorus, the values of the flux are actually comparable with the values for 
the KAM for small periods, but the last we were able to compute is of the order of 10~ 9 . 



The Fig. (24) shows that it seems to already converge towards a value larger than the one 



obtained for the KAM. 

In conclusion, the Greene's conjecture and Mather's theorem are verified in both situ- 
ations and the numerical results are in agreement with the theory. From Fig.(p4|) we can 
thus conclude that the existence of (i) a KAM surface, the plasma edge barrier, having the 
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Flux 




m 

Figure 24: Comparison of the decrease of the flux across convergent island chains towards 
the upper Cantorus q = N(l, 7) (circles) and towards the KAM surface on the plasma edge 
q = TV (4, 2) (squares). The former is seen to converge to a finite value of the order of 
10 -8 , while the latter is observed with lower values to continue to decrease hopefully to zero 
(robust barrier). 
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rotation number N ^ 2 ) anc ^ °^ fa) a P^tial barrier, a Cantorus, with the rotation number 
N ^ 7 ^ was indeed confirmed by using numerical methods. 

6 Asymptotic radial subdiffusion 

Dispersive motion of magnetic lines in a stochastic zone is often described very roughly as 
a random walk. We stress the fact that the radial motion observed here is non- diffusive. 
In order to determine more precisely the average dispersion properties of magnetic lines in 
a situation of incomplete chaos (diffusion, subdiffusion or superdiffusion) and their scaling 
properties with the perturbation, we have calculated the asymptotic properties of a set of 
5000 magnetic lines. Initial conditions are chosen in the chaotic shell on an initial circle 
ip = 0.001, r ~ 0.0316, outside of the protected plasma core, either as a bunch of lines 
within a small angular interval A9 = 0.001, or on a complete circle with regularly spaced 
initial values of the initial angle 9 between and 1. The values of the stochasticity parameter 
L are varied in a wide range between 0.776 and 6.68, beyond the threshold for large scale 
chaos above which all magnetic lines from the chaotic shell finally cross the edge barrier, 
transformed into a permeable Cantorus. Here we also study the radial position x = yfip of 
a magnetic line, which is unity on the unperturbed plasma edge. 

We have to stress the fact that with such large values of the stochasticity parameter 
L, almost all magnetic lines in the central chaotic shell are crossing the plasma edge. This 
situation is of course of less interest for confinement studies. It would rather be characteristic, 
instead, of a situation of plasma disruption due to the rapid escape of the magnetic lines. The 
time behavior of this escape is however interesting in order to study the scaling properties of 
this dispersive motion, as compared with other maps like the standard map, where a similar 
flux diffusion has been intensely studied since many years. 

6.1 Long time behavior 

Although the individual motion of one magnetic line appears to be discrete and random, 
average properties like the mean poloidal flux tp m and the mean radius, averaged over initial 
conditions, can be described by "continuous" functions of "time" (the number of large turns 
around the torus) : 

1> m (t) = ty{t)) , x m (t) = (\/#)) , (94) 

In order to deduce dispersion properties, we also calculate a third independent quantity : 
the mean square displacements (MSD) of the flux, defined in the reference frame of the 
average motion : 

MSD^(t)^(m)-An(t)} 2 ) (95) 

The MSD of the radius is also computed, and we check that the following exact relation is 
verified at each time in terms of the average motion : 

MSD x {t) = ^ [VV^) - %(*)] J = iM*) - xm(t) 2 (96) 

The intuitive picture emerging from these results is the dispersion of magnetic lines from 
an initial circle of intersection points, passing local weak barriers, coming back an forth and 
expanding radially in average. The asymptotic state of course strongly depends on the value 
of the stochasticity parameter L , as well as all characteristic times which appear to scale 
like L~ 2 (see Eqs.fllO^) below). 
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6.1.1 Asymptotic motion in the confined domain 



In the confinement domain, the global diffusion of course vanishes, but the details of the time 
behavior of an ensemble of magnetic lines is characteristic of the phase portrait presented 
in Fig. dTT|) . For L = 4.875/2-7T ~ 0.776, the orders of magnitude of the characteristic times 
can be inferred from Fig.(|l^) and tested in the present average evolution : residence times 
inside the ITB would be of the order of tjtb ~ 10 6 (downwards peaks), residence times in 
the chaotic shell would be of the order of Tg^eii w 10 7 and in the chaotic sea of the order of 
T sea ~ 10 8 iterations. 



r itb ~ 


10 6 


T shell ~ 


10 7 


'Tsea ~ 


10 8 



(97) 



We first consider a bunch of 500 initial conditions within A9 = 0.001 followed over 10 9 
iterations for L = 4.875/2-7T ~ 0.776. We observe that none of the 500 considered magnetic 
lines succeeds to escape from the plasma edge even after 10 9 iterations, and the average 
radius saturates around a value x m (t) ~ 0.7 which proves that the robust KAM barrier on 
the plasma edge is still resistant (see Figs.(fl3|) and ([ll])). The detailed form of the time 
dependence is however interesting, especially concerning the average radius x rn (t) and the 
radial dispersion MSD x {t) (see Fig.(|2^). In Fig.(p5|.a) which represents x m (t), we see at 
short times large oscillations starting from the initial condition x ~ 0.0316 and saturating 
after about 10 2 iterations; this can be explained by the motion of the bunch of lines as a 
whole, along the non circular magnetic surfaces making an excursion towards larger x values 
at each poloidal turn in 9. 

For larger times, up to 10 6 , a slow increase of the average value is observed up to a value 
x m (t) ~ 0.2, while the dispersion remains constant (Fig.(^.b)). This value x m (t) ~ 0.2 is 
important since it roughly corresponds to the position of the transport barrier (averaged 
along 8). The resulting interpretation is the following : the ensemble of magnetic lines 
describe a CTRW with long sojourn times inside (or below) the ITB, so that the average 
position does not exceed that of the ITB. The dispersion remains constant (zero diffusion) 
since the chaotic motion occurs inside a restricted region below the upper Cantorus of the 
barrier. 

Understanding the global motion for longer times, between 10 6 and 10 8 , appears as 
a challenging question. It is seen indeed in Fig.(|25|.a) that the average radius suddenly 
begins to increase in time (roughly like t 1 ^ 2 ), while the dispersion MSD x (t) (Fig.f^.b)) 
first increases (roughly like t 1 ) from 10 6 to 10 7 , has a maximum and then decreases toward 
a constant value. An interpretation of this overshooting can be tentatively proposed as 
follows. 

- For times shorter than titb ~ 10 6 the dispersion is dominated by the trapping processes 
inside the ITB, and the very few lines which succeed to escape across the barrier actually 
create a slow increase of the curve. 

- Then from tjtb ~ 10 6 to r shell ~ 10 7 usual diffusion occurs in the chaotic shell and 
MSD x (t) grows linearly in time. 

- Around T s h & ii ~ 10 7 saturation of the dispersion is observed because of the finite size 
of this chaotic shell limited by the ITB. 

- After T s heii ~ 10 7 a non negligible part of the magnetic lines are in the chaotic sea, 
and long sticking events occur around the main rational chains (see Fig.(p2|)) which result 
in a decrease of the dispersion. 

- Final saturation of the dispersion is obtained around r sea w 10 s with an average value 
x m (t) -0.7. 
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AVERAGED NORMALIZED RADIUS x DISPERSION OF THE NORMALIZED RADIUS 




Figure 25: Average radius (a), mean square radial dispersion (b), average flux (c) and mean 
square flux dispersion (d) of an initial bunch of 500 lines at L = 4.875/2-7T followed during 
10 9 iterations. 
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This average position of the magnetic lines can also be roughly explained in terms of 
the estimated average residence times : each region (chaotic shell, ITB and chaotic sea) has 
indeed an average radial position (averaged over 9), and when these average positions are 
weighted by the percentage of time spent in each region, one finds indeed a result of the 
order of the asymptotic radius observed in Fig.([25[a): x rn (t) — > 0.7. Such an agreement 
should not be taken too seriously, and very long iterations should be necessary in order to 
compute accurate residence times for this presumed CTRW, but the orders of magnitude 
found yield a first explanation. 

When the same calculation is performed for an initial circle of magnetic lines starting 
from the same value ip = 0.001, r ~ 0.0316, analogous results are obtained in Fig.(|26|) except 
for the short time oscillations which are of smaller amplitude but which however remain. 



6.1.2 Asymptotic motion in the transition domain 



For a slightly higher value of the stochasticity parameter L = 5/2tt ~ 0.796, the same 
series of transitions occurs at small value of the "time" : we may thus consider shorter 
asymptotic times and a larger number of initial conditions. For such L value, all magnetic 
lines actually go beyond the plasma edge, indicating that the KAM on the edge is already 
broken, transformed in some Cantorus. The latter still represents some barrier, and the 
average motion is different inside and outside this permeable barrier. 



In Figs. (||^. a to d) we have represented the four functions x m (t), MSD x (t), ip m (t) and 
MSD^(t) for a bunch of 5000 initial lines. Clearly a similar transition region is obtained, 
as in Figs. ( p5| ) and (p6~|), but here for smaller values of time (roughly from 10 5 to 10 6 ). In 
this time domain a straight line is observed in each graph, corresponding to the following 
transient time behaviors: x m (t) ~ i 1/2 , MSD x (t) ~ t 1 , tp m (t) ~ t 1 and MSD^(t) ~ t 2 
which represent a transient diffusion of the radius and a superdiffusion of the flux ip. For such 
an initial bunch of lines, the initial oscillations are well pronounced. Here the asymptotic 
behavior is not yet reached after 10 8 iterations and it seems that the dispersion of the radius 
MSD x {t) could tend to a constant value, probably because of the presence of a resistant 
barrier very far out of the plasma. The latter is of no physical relevance for the tokamak, 
but is however characteristic of the present process of destruction of successive barriers in 
an incomplete chaos situation. 



6.1.3 Asymptotic motion in the large scale dispersion domain 

Beyond such transition region, for L ~ 5.5/2n ~ 0.875 the time behavior is more regular. In 
Fig. (28. a to d), we have represented the four function x m (t), MSD x (t), ip m (t) and MSD^(t). 
As expected the initial oscillations are less pronounced, but still exist. Some transient 
behavior can hardly be seen (around 10 3 in MSD^(t)) and a rapid growth leads to a final 
asymptotic state reached after 10 5 — 10 6 iterations, in which a simple time behavior is 
measured. 



One notes a common (but not fully understood) feature of all the graphs in that domain 
of i-values: all these functions seem to present an inflection point when reaching the absolute 
value unity. In the graphs for x m (t) this occurs at a time of physical relevance, called the 
escape time t v (the time to reach the edge of the plasma at tj) = 1, x = 1, such that 
x m (t v ) — 1) which is here t v ss 10 4 . 
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AVERAGED NORMALIZED RADIUS x DISPERSION OF THE NORMALIZED RADIUS 




Figure 27: Average radius (a), mean square radial dispersion (b), average flux (c) and mean 
square flux dispersion (d) of an initial circle of 5000 lines at L = 5/2ir. 
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In Fig. (^8|.b) for M SD X (t) we remark that some overshooting of the curve appears around 
10 5 , before the asymptotic state is reached. Such bump does not appear on the other func- 
tions. Is it simply the effect of a distant barrier, out of the plasma edge, which temporarily 
traps magnetic lines ? 

Such a bump also appears in Fig.(p9|.b) for MSD x {t) at a higher value L = 6/2ir ~ 0.955 
for an initial bunch of 5000 lines. In each of the Figs.(|2^.a to d), all the features remain 
unchanged, except that the escape time is decreased i„ f=a 10 3 and that asymptotic state 
is reached more rapidly, around 10 . 



For higher L values L — 9/27T, dispersion is still more rapid and the asymptotic state 
is reached more rapidly. Here we follow a circle of 5000 initial lines. One note however in 
Figs. (30 , a to d) that the curve in log-log scale has a long time slope which is slowly decreasing 
towards what we consider as the asymptotic state, so that long iterations are necessary to 
evaluate accurately the time behavior and the dispersion coefficients introduced below (see 
Eqs. (|02|) ). 



For a set of values of the stochasticity parameter L in the large scale dispersion domain, 
we have measured the asymptotic slope and dispersion coefficient of each quantity (^ to 
^6|) . Our measurements indicate that the average radius is found to increase asymptotically 
in time as t 1 ^ 4 : 

x m (t) -> a t 1 ' 4 (98) 

In a consistent way we find 

^ m (i) 6 tV* , MSD^t) -> 2 A/, t (99) 

i.e. diffusion of the flux, like in the standard map, and of course 

MSD x (t) —>2D X i 1/2 (100) 

which means radial sub-diffusion. Although surprising, this is however not a fully new result: 
similar phenomena of "strange" transport 14 p3] (where MSD r (t) ~ f with fx =/= 1) have 
been measured, namely for charged particles sticking around an island remnant in a chaotic 
three- wave chaotic model |HJ resulting in a t 1 / 3 behavior p2| . 



It has to be noted that, with the same unperturbed magnetic field (p5|), the direct 
integration of the magnetic line equations of motion, in presence of Fourier series magnetic 
perturbations or a distribution of current filaments, has already yielded a fully equivalent 
i 1 / 4 asymptotic behavior the quantity y/MSD x (t) rapidly called by these authors as being 
some "average displacement " of magnetic lines in chaotic layers This W 4 law is in total 
agreement with the results ( |100[ ) obtained here in a simple map. 



The above subdiffusive result (100) is however reached asymptotically for long times. In 
the present incomplete chaotic regime of the tokamap, it shows a slower dispersion (fj, < 
1) which is different from the usual diff usio n in a completely chaotic situation [pi[ . A 
comparison is performed below, after Eq.(112) and in Fig.(|32]), in Section 3.4 . 



For physical times, of the order of the average time t v necessary to escape from the 
boundary circle, different transient regimes are observed instead, with a rapid radial super- 
fusion (/j > 1) observed when crossing the plasma edge. In other words, although the 
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AVERAGED NORMALIZED RADIUS x DISPERSION OF THE NORMALIZED RADIUS 




AVERAGED FLUX versus TIME ( LogLog ) ==> DISPERSION OF THE FLUX versus TIME ( LogLog ) 
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Figure 29: Average radius (a), mean square radial dispersion (b), average flux (c) and mean 
square flux dispersion (d) of an initial bunch of 5000 lines at L = 6/2n. 
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AVERAGED NORMALIZED RADIUS x DISPERSION OF THE NORMALIZED RADIUS 




AVERAGED FLUX versus TIME ( LogLog ) ==> DISPERSION OF THE FLUX versus TIME ( LogLog ) 
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Figure 30: Average radius (a), mean square radial dispersion (b), average flux (c) and mean 
square flux dispersion (d) of an initial circle of 5000 lines at L = 9/2ir. 
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asymptotic, long time properties are diffusion of the flux and sub-diffusion of the radius, 
several other behaviors are observed along the time evolution, with a slower or faster time 
variation. 

The numerical results measured for the coefficients a, 6, D x and are the following, 
as function of the stochasticity parameter L : 
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(101) 



6.2 Scaling properties 

Another important feature of transport processes is the scaling property of the corresponding 



dispersion coefficients (101) in terms of the perturbation parameter L. In order to draw the 
scaling laws for these coefficients, we define arbitrarily the first values by L a , a , b , D xo and 
D and we plot the above numerical results reduced to these first values as L/ L Q , a/a a , b/b , D x /D xo 
and D^/Dtpq, respectively. 

We clearly see in Fig.(|l|) that the various coefficients scale as follows : a(L) ~ L 1 ' 2 , 
b{L) ~ L ~ D X (L) and D${L) ~ L 2 . 

Summarizing the results obtained for the both dependences in L and t (see Section (@)), 
we can write down the following asymptotic formulas describing both the long time depen- 



dence (Eq.(|98|) to (100)) and the scaling properties of the so-defined dispersion coefficients 



a(L), b(L), D X (L) and D^L) 



1/4 



lim^ x m (t) = a(L).t 1 /4 ~ [£2 t ] 

lim^co iiJ m (t) = 6(£).tVa ~ [ L H] 1/2 (1Q2) 

lim^oo MSD x (t) = 2 D x (L).t 1 / 2 ~ [lH^' 
lim^ MSD^t) = 2 D^(L).t ~ [L 2 i 

This clearly shows a global scaling in \_L 2 t\ . Such scaling laws can be simply interpreted as 
resulting from a simple intuitive law of the type 5if)(L, t) fa [fix 2 ] — L. 

In the standard map, which is periodic in both variables angle and action, it is well 
known ](||, Q (see also p. 299 in Q) that computation of the diffusion coefficient of the 
action yields the result D(L) = D q i(L).S(L), where S(L) involves a series of Bessel function 
J(£), and where the simple coefficient 

D ql (L) = L 2 /A (103) 

can be obtained from a random phase argument. This classical result for the standard map 
is actually modified by the presence of accelerator modes, which do not occur here. Beyond 
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Figure 31: Graphical presentation of the scaling dependence in the stochasticity parameter 
L of the asymptotic coefficients of flux diffusion ~ L 2 (black circles) , of radius subdif- 
fusion D x r~j L (grey squares), of average flux growth b ~ L (black triangles) and average 
radius growth a ~ L 1 / 2 (black diamonds). For a lowest value Lq — 5.25/2-7T ~ 0.836 this 
graph represents these various coefficients C in a logarithmic plot of log(C/Cb) as function 
of log(L/Lo). The three straight lines correspond, from top to bottom, to the expected 
behaviors in L 2 , L and L 1 / 2 , respectively. 
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the scaling relation D^(L) ~ L 2 (Eq. ( |102| )) found to be satisfied in the tokamap, it is 
interesting to compare the exact numerical results with Eq. (|103|). We find 



D^(L) ~D ql (L). [1±0.1] 



(104) 



which shows not only that the same quasi-linear diffusion coefficient D q i(L) applies both 
in the standard map and in the tokamap, but moreover that no Bessel functions seem to 
appear here, probably due to the lack of periodicity of the map in the flux variable. 



Some methods for the calculation of diffusion coefficients in maps where variables are 
clearly separated are given in p7j. An analytical derivation of even a simple result like 
Eq.(104) for the tokamap (Eqs. (|31|), ( |32| ) where variables are not separated) seems however 
to be difficult to prove in general using the usual Fourier path diagram method 66 . It is 



however very simple to check that the same random phase argument can be applied to the 



tokamap and the same result (103) is obtained, as for the standard map but without any 
Bessel functions here. Unlike the standard map, this asymptotic behavior (104) is valid for 
all values of the stochasticity parameter L. 



6.3 Asymptotic times and transient regime inside the plasma be- 
yond escape threshold. 

The asymptotic state for L = 5.5/27T ~ 0.796 (beyond escape threshold) has been obtained 
after a rather long relaxation times 

t R » 10 6 . (105) 

We have no theoretical explanation for such slow phenomenon in the tokamap, but in similar, 
simpler systems, relaxation times are also very long. For the standard map j2^] p3[ , very 
long iteration sequences have been necessary to perform measurements of the transition 
probability between different basins and to characterize the corresponding CTRW. A detailed 
kinetic theory has been elaborated very recently |58j for the standard map, showing also 
very long relaxation times. 

Since asymptotic behavior is only reached after such long relaxation times, of course 
magnetic lines are already far beyond the physical domain of the plasma (tp ^ 1). For 
instance in Fig.(|2g|.a), for a value L = 5.5/2n just after breaking of the edge KAM, we 
can see that the average radius r m reaches the edge of the (unperturbed) plasma after a 
"physical" or "escape time" t v such that r m (t v ) = 1, which yields 

t v « 10 4 (106) 

At least for values of L > 5/2ir ~ 0.796 , beyond the escape threshold, this time represents 
an estimate of the escape time, or more precisely, an estimation of the number of turns (the 
long way) around the torus after which a magnetic line has reached the plasma edge. 

It is interesting to translate the number of iterations in the tokamap (thus the number 
of turns) into an average time for thermal particles to reach the plasma edge by following 
the magnetic lines (in absence of collisional or magnetic drift effects). One iteration time 
step in the tokamap represents actually the length of one large turn of lengt h 2ttR q around 
the torus. For a thermal particle following the field line the physical time (|l06| ) can thus 
be evaluated as the time necessary for a particle to travel t v times the long way along the 
torus, i.e. 2ttRq * t v /Vth- In Tore Supra, this time is of the order of 10~ 2 s. for ions and 
10 -4 s. for electrons in usual conditions (T = 5 keV) 

t vi = 10- 2 s. , ^ e = 10- 4 s. (107) 
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In the same way one can deduce that the relaxation time, corresponding to 10 6 iterations 
steps of the tokamap, actually represents 20 s for ions and 0.5 s. for electrons in Tore Supra. 

t R i = 20 s. , t R e = 0.5 s. 

These are the times necessary to reach the slow asymptotic subdiffusion regime at L = 
5.5/2-7T ~ 0.796 . The transient superdiffusive regime occurring for t < t R is thus not 
excluded from the domain of practical interest, as well as asymptotic subdiffusion regime 
obtained here. 



These values of the escape time in the tokamap (107) are thus much shorter than the 
'particle confinement time Tp >z te ~ 0.20 s. which is admitted to be roughly equal to the 
energy confinement time te, m the absence of additional heating for Tore Supra for instance. 



Such escape times ( 106 ) are thus relevant for particle confinement losses in tokamaks. 



For times of the order of the escape time t j§ t v> the average radius is growing like 
x m (t v ) ~ t 1 ', much faster than in the asymptotic state (p8|). The general shape of the curve 



Fig. (28. a) indeed indicates a slow initial growth, followed by a rapid increase of the slope, 
similar to an exponential growth, and a final saturation with a lower asymptotic slope. We 
note that the pseudo-exponential growth occurs precisely around the physical time (|106|). 



In other words, although the asymptotic behavior of magnetic line motion in the tokamap 



is a rather regular one (even if a "strange" one) described by (|102| ), nevertheless for finite 
times of the order of the physical time we observe the following extremely rapid transient 
time-behavior : 

MSD^{t 9 ) -t 5 - 88 

MSD x {t v ) ~t 2 75 (108) 

3>m(fip) ^ t 

In the transient regime occurring inside the plasma, for L — 5.5/2-7T, we have thus a su- 
perdiffusion MSD x (t v ) ~ i 2 - 75 of the magnetic lines, around an average radius growing like 
x m (t v ) ~ t i.e. in a ballistic way. This last observation would seem a priori much less 
favorable for plasma confinement in perturbed magnetic structures, but it is actually not 
the case. 



6.4 Comparison with classical diffusive predictions 



We have to compare the present result (108) with the classical analytical prediction, i.e. 



the diffusion coefficient of magnetic lines in a completely stochastic magnetic field. This 
process is parametrized by the amplitude of the magnetic perturbation [3 =| SB \ / B and by 
parallel and perpendicular coherence lengths of the perturbation, A M and Aj_, respectively. 
This naturally introduces the control parameter of the problem : the Kubo number 

a = P^- (109) 

which roughly measures the parallel coherence of the fluctuating magnetic field : magnetic 
diffusion will be small if a is small. The quasilinear approximation Dql for the diffusion 
coefficient in a completely chaotic field is simply : 



DQL = T^-\lo- a2 = \hP 2 K (no) 



It is actually a quadratic approximation in the Kubo number a , where the first factor 
takes the dimensionality into account (the units of the magnetic diffusion coefficients are in 
m 2 /m.). Higher order terms have been derived from a more general description of diffusion 
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in a completely stochastic fields in [Q. This coefficient (11C) has been originally derived 
in p9[ , J70| . The radial mean square displacement < x 2 (z) >= 2Dql.z is growing linearly 
with the toroidal distance z . 

In order to make a connection with physical parameters of the tokamak, we have to 
know the relation between the stochasticity parameter L and the relative strength of 
the perturbation introduced for instance in the non-axisymmetric external coils. From the 



alternative derivation of the tokamap from canonical variables 32 1, this relation can be 
derived as : 

-1(f) 

where Ro/a = 1/st is the aspect ratio of the torus. 

Within the present dimensionless units (reduced to the minor radius a of the plasma 
column), and with t being as previously the number of long turns around the torus, this 
classical result writes (with z = 2TrRo.t, A, » 2ttRq, r 2 = 2x 2 and L = 2-Kjjjet) : 

MSD QL (t) = V2^L 2 t (112) 



It can easily be shown that this linear function ( |l 12| ) remains for all times much above the 
numerical result of Fig.(|2q.b), even for physical times of the order of t v where the growth 
rate is like t 2 75 : in Fig. (p2[) we draw the dispersion measurements of the radial position in 
the present situation of incomplete chaos for L = 5.5/27T ~ 0.876 (as in Fig. fl2S|.b)) along 
with its asymptotic slope (dotted line) which exhibits an asymptotic radial subdiffusion. In 
the domain between t = 10 3 and 10 4 (around the escape time t v , see Eq. ( |106[ )), one clearly 
observe a superdiffusive regime. All these measured behaviors of MSD x (t) nevertheless ap- 
pear at all relevant times to remain much smaller than the quasi-linear diffusion in complete 



chaos as represented by the upper straight line described by the classical Eq. (112) |6 



In conclusion, the transient superdiffusion observed in the tokamap model for times of 
the order of the physical or escape time is actually less effective than the classical quasilinear 
prediction. This means that classical transport in a completely stochastic magnetic field as 
studied in |^|, |70| and is actually more dramatic for plasma confinement than the 
strange transport observed in the present model for incomplete chaos. A transient, very 
rapid growth, although "superdiffusive" is not always "more dangerous" than a linear one... 



7 CONCLUSIONS: 

The chaotic motion of magnetic lines in a toroidal perturbed geometry has been analyzed 
using a Hamiltonian map called tokamap [jl6| defined by (|32|, [33[ |34|) in terms of the stochas- 
ticity parameter L. It describes the successive intersection points of a magnetic line in a 
poloidal plane (ip,6) . From the previous analysis of fixed points bifurcations, new results 
have been derived for the Shafranov shift in this map (Fig.(|l|)). We restrict ourselves to the 
case of an monotonous unperturbed safety factor profile with a central value q(0) = 1. 

The generic phase portrait in the poloidal plane (tp, 9) consists in a protected non- 
chaotic plasma core, surrounded by a chaotic shell limited by a semi-permeable Cantorus 
which has been identified. For increasing values of L magnetic lines finally cross this first 
barrier and wander in a wide chaotic sea, extending up to the plasma edge (Fig.(]TT|)). This 
domain of the phase portrait has the realistic structure with remnants of island chains on 
the dominant rational q- values (Fig.(^j)), and sticking process of magnetic lines around these 
island remnants. 

The important point we want to stress is that the different barriers found in this con- 
figuration always correspond to irrational values of the safety factor, and not to rational 
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Figure 32: Dispersion measurements of the radial position in the present situation of incom- 
plete chaos are represented (lower black curve) for L = 5.5/27T ~ 0.876 (as in Fig.Dif 3b) 
with its asymptotic slope (dotted line) which exhibits an asymptotic radial subdiffusion. In 
the domain between t = 10 3 and 10 4 (around the escape time t v , see Eq. 99), one clearly 
observe a superdiffusive regime. All these measured behaviors of MSD x (t) nevertheless ap- 
pear at all relevant times to remain much smaller than the quasi-linear diffusion in complete 
chaos as represented by the upper straight line described by the classical eq. 105 of Ref.[69] 
(Rechester & Rosenbluth 1978). 
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values where island chains are formed. The golden number G is known to play no role in 
this map p7[ |. It is thus quite natural from KAM theory that the "next most irrational" 
numbers, the "noble" numbers N(l, k) = 1 + l/(fc + G" 1 ) jfTJ could actually play an im- 
portant role. The wandering motion in the inner shell is indeed limited by a KAM surface 
(or a robust Cantorus) which is located in the case L = 4.875/27T ~ 0.776 downward at 
q = N(l, 11) = 1.086.., and upward by a semi-permeable Cantorus at q — N(l, 8) = 1.116.... 
The wandering motion in chaotic sea is limited downward by another identified Cantorus, 
which is located slightly above the previous one at q = N(l,7) = 1.131 : we thus find the 
existence of a double-sided internal transport barrier. Although this barrier includes a ratio- 
nal surface (q = 9/8 = 1.125 in the case L — 4.875/27T ~ 0.776) we have shown that the two 
Cantori, which are actually the two sides of the barrier, are located on noble values of the 
safety factor. The outer edge of the chaotic sea remains a strong external transport barrier 
(a KAM surface or a robust Cantorus) at L = 4.875/27T ; its radial position in the q -profile 
appears to be located with a rather good precision at the most irrational number between 
q = 4 and q = 4.5, i.e. at q = A(4,2) = 3 + A(l,2) = 4+l/(2 + G~ 1 ) = 4.382... The 
transport barriers observed in this map have been proved (from flux argument) to be noble 
Cantori. The flux of lines across permeable Cantori has been discussed in Section ^, by 
using very precise calculations of high order periodic points along rational chains convergent 
toward the Cantorus. The method used for determining these periodic point is summarized 
in Appendix B. 

The motion along a very long trajectory shows moreover that the magnetic line motion 
is intermittent ( Fig.([13|)), with long stages in the inner shell and in the outer chaotic 



sea, with occasional sticking periods around island remnants (Fig. (12)). These features are 
characteristic of a continuous time random walk (CTRW). The external barrier actually 
confines magnetic lines inside the plasma up to the appearance of a large scale motion above 
a threshold of the stochasticity parameter around L ~ 5/2-7T. 

In order to describe the motion of magnetic lines beyond this large scale motion threshold, 
we have considered a bunch of magnetic lines in a small area of phase space and study their 
statistical average evolution. We have computed, the average radial position (|94|) of the 
bunch of lines, and the mean square deviation of the radius ( |96| ) compared to this average 
radius : this is the quantity defining the spatial diffusion properties in the asymptotic stage 
(subdiffusion or superdiffusion) . We have also calculated the diffusion ( p5| ) of the poloidal 
flux ip : the latter appears to describe asymptotically a classical diffusion (in if> !) and behaves 
as traditionally admitted, i.e. like the standard map in the quasi-linear stage, with a diffusion 
coefficient in action space which scales like D ~ L 2 in the domain L — [5/2-7T, 21/2vr]. This 
diffusive behavior does not hold for magnetic lines in physical space. 

Since the action of poloidal flux ip is proportional to the square of the radius x, this 
diffusion process in action is equivalent to a fourth moment of r. It is thus not surprising 
at all that the radial r-motion appears to be asymptotically a spatial subdiffusion 

< \x(t) - x m (t)} 2 >=> 2D X (L) t 1 ' 2 with D X {L) ~ L (113) 

In contradistinction with usual results for dispersion in completely chaotic situations (quasi- 
linear description) , the radial motion of stochastic magnetic lines in the present tokamap 
model, with incomplete chaos, thus appears as a radial "strange diffusion" and not classical 



diffusion. In (113) x m {f) is the average radius 



x m (t) =< ^/W) >=^ a(L) t 1 / 4 with a{L) ~ L 1 ' 2 (114) 
which grows asymptotically as (L 2 t) 1 / 4 . The scaling laws (102) for the coefficients of these 



three average processes reveal a nice scaling in (L 2 t) for the asymptotic diffusion properties 
of the tokamap, even for strange diffusion. 

These general and asymptotic properties of magnetic line motion in a perturbed toroidal 
geometry have been obtained in the tokamap model. However, the threshold value of the 
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stochasticity parameter L *~ 1 , for which particles begin to escape from the plasma to 
describe large scale motion and a true diffusion in action space, is sufficiently high, and the 



relaxation time tji (105), necessary to reach the asymptotic stage, is sufficiently long for this 
asymptotic stage to be reached when the average radial position of the bunch of lines has 
already crossed the plasma edge ip = f . As a consequence, the magnetic line motion inside 
the plasma is either confined by a KAM external barrier, and the asymptotic motion can 
only be subdiffusive when L < 1 (like in Ref. p3[), or escaping from the plasma when L > 1 , 



but this occurs in average at a physical time ttp (|106|) much shorter than the relax ation time 



tff , so that the motion inside the plasma is not yet the asymptotic one (see Eq. 108), and 
actually describes a transient regime of radial superdiffusion. This last observation could 
seem much less favorable for plasma confinement in perturbed magnetic structures with 
incomplete chaos when L > 1, but actually the resulting motion appears to be much less 
dramatic for plasma confinement than the usual classical expectation from the quasilinear 
description of classical diffusion in completely stochastic magnetic fields [f70| , |p4|| . 
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8 APPENDIX A: Standard safety factor profile. 

Standard radial profiles can be derived, within a circular cross section, by taking into account 
elementary observations in cylindrical geometry (p, 9, z). The basic assumption of the model 
consists in assuming a parabolic density profile with vanishing condition on the edge: 

n(x) = n(0) [l - x 2 ] (115) 

where the radius x = p/a is reduced to the small radius a of the plasma. Such a profile is 
not inconsistent with early observations of ohmically heated plasmas fnj . 

Next, a general relation is assumed between temperature and density profiles. A system- 
atic fitting of many profiles on various early tokamaks ]72|] indicates that the temperature 
profile can adequately be represented by squaring the density profile. Such a result can also 
be obtained from simple analytical models of energy balance (see e.g. p. 41 in (73)). We 
thus consider : 

T{r) =T(0) [l-a; 2 ] 2 (116) 

which also implies r\ e = d In T e /d In n e = 2. 

In order to obtain the density current profile, we assume a Spitzer dependence of the 
resistivity i] ~ T 3 / 2 as function of the temperature j74| . For a stationary magnetic field, 
Faraday's law implies indeed an irrotational electric field ^ = -VxE = hence in an 
axisymmetrical system = 0, i.e. a constant electric field along the radius. From a 
simple expression of the Ohm's law we thus find that E z (p) = rj{p)j z {p) = est. is actually 
a constant along the radius p. From the temperature dependence of the Spitzer resistivity 
one thus deduces in this simple model that the electric current density profile is 

j z (x)=j z (0)[l-x 2 } 3 (117) 



where z is the toroidal direction in cylindrical geometry. We note that Eqs. (115, |116| and 

(118) 



117) imply 

{3z)j, (p), 



jz(0) p(0) 
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where p is the plasma pressure and (...) s denotes the average inside a magnetic surface tp- 



This relation (118) has been demonstrated to hold in ohmically heated discharges of the 



TCV tokamak, as quoted by Minardi |75| . He proved that the relation (118) follows from 
the stationarity of the magnetic entropy, in general conditions where the ohmic tokamak is 
considered as a "dissipative open system in equilibrium which interacts externally with the 
ohmic transformer" . 

We may then use Ampere's law in cylindrical geometry 

j z = (VxB) z = ~(pB e ) (119) 



and integrate the current density profile (117) to obtain 

Bg{x) = B e {l)x (2 - x 2 ) (x 4 - 2x 2 + 2) (120) 

where the value Bg(l) on the edge is given by the total electric current / in the plasma : 
B e (l) = aj z (0)/8 = I/2ira. 

In the standard model for the toroidal magnetic field the safety factor q(p) is defined 
in terms of Bg and on the central value -Bo of the toroidal magnetic field on the magnetic 
axis by q(x) — etxBq / 'Bg(x) so that we finally obtain the safety factor profile in cylindrical 
geometry as 

^ = (2-, 2 )(ji 0) 2, 2 + 2) ^ 

This standard profiles indicates a value on the edge four times larger than on the axis 
q(l) — iq(0) which is rather reasonable in most usual ohmic cases. By denoting x 2 =>- ip, we 
obtain the expression actually used in |]l6| : 

m = pP (122) 

K ' (2- V)(^ 2 -2^ + 2) V ' 

where -0 varies from in the center to 1 on the edge. 

In order to make use of traditional canonically conjugated coordinates (r, #), see Eq. (|36|) , 
we may also introduce r 2 /2 = x 2 

^ = - (123) 

where the modified coordinate r = x\[2 varies from in the center to \/2 on the edge. The 
difference between x = p/a and r = p\/2/a may simply be interpreted in terms of the r as 
occurring in a plasma of radius \[2. In either way the edge of the plasma is located at tp = 1, 
x = 1, r = V2 and we generally use plots in terms of the radius x varying from to 1 inside 
the unperturbed plasma edge. 



9 APPENDIX B: Numerical method for finding peri- 
odic points in discrete maps 

In this appendix we explain how to determine hyperbolic and elliptic periodic points, by a 
numerical algorithm derived from a generalization of the Fletcher-Reeves method, involving 
the Jacobi matrix of the tokamap. More details are given in |f76f . 
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The existence of points with everywhere dense trajectories and the possibility to ap- 
proximate every point by a periodic point are two important and apparently contradicting 
aspects of chaotic systems. Nevertheless, the localization of many periodic points appears 
highly necessary in order to describe the structure of the phase space. This is because a typ- 
ical Hamiltonian system appears to belong to an intermediate category between completely 
integrable and globally chaotic systems. Their phase space decomposes in a well organized 
architecture of invariant subdomains. In this architecture the periodic points play the role 
of a skeleton. The knowledge of their relative positions and residues, permits to determine, 
at least qualitatively, the whole phase space structure. 



The residue's absolute value (see Eq.(76)) of hyperbolic periodic points increases expo- 
nentially with the order of periodicity m. This exponential increase imposes clear restrictions 
on the practical possibility of numerical localization of the higher order hyperbolic periodic 
points. The localization of periodic points approaching an invariant circle is facilitated by 
a - let us say a continuity principle, as follows. Consider the residue of a sequence of chain 
of elliptic or hyperbolic points, convergents of an irrational, generic, KAM barrier : then 
the residues of periodic orbits will ap pro ach the residues of the invariant circle, i. e. will 
approach to zero (Greene's conjecture |L4j ). In any neighborhood of an invariant circle there 
are elliptic and direct hyperbolic points. 

A quite inverse phenomenon occurs at half-pcrmcable barriers, the Cantori, where the 
typical values of the residues tend to -co for direct hyperbolic points and to +00 for inverse 
hyperbolic points. 



9.1 Methods for computing periodic points (theoretical aspects) 

We precise that in the following, the angular coordinate (denoted by x\) will not be a 
priori reduced modulo 1, despite the fact that in the numerical calculations, for an increased 
numerical accuracy, we separately represent the integer and the factional part of the angular 
variables. After each iteration, the integer and fractional parts of the angle are computed 
and stored separately. 

For our purpose, in order to find a periodic point of order of periodicity m and q = m/n, 
we must solve a system of a nonlinear equation, like 

Xl = TM 1 (xi,x 2 ) -n (124) 
x 2 = TM 2 { xt,x 2 ) 

where TM = T m is the m times iterate of the original map T. The first equation simply 
expresses the fact that the m th iterate TMt{x\, x 2 ) of the angular coordinate has the same 
value as the initial one (xi) but increased by the integer number n of periodic points. 

The problem of localizing periodic orbit is reduced to solving a system of nonlinear equa- 
tions. 

More generally, we can consider the following problem: elaborate an algorithm for solving 
the system of N nonlinear equations 

F i (x)=F i (x 1 ,x 2 ,..x N ) = (l<i<N) (125) 

by using the values of the functions Ft, F 2 , F^. 

There are a lot of numerical methods for solving systems of nonlinear equations. Because 
of the large roundoff errors we are obliged to choose a very rapidly convergent method. 
We chose the Fletcher-Reeves minimizing method because it has a very rapid convergence 
(double exponential) and a basin of attraction which is large enough. 
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Well known iterative methods seem to be not appropriate for our aim. Consider for 
instance the traditional Newton iterative method fin\ . It is known to converge quadratically 
(near the exact root, the number of significant digits approximately doubles at each step) if 
the initial point is chosen in the basin of attraction of the root, but the global convergence 
properties are very poor and unpredictable [ [T7| . The basin of attraction is usually a fat 
fractal set whose boundary can be very near the root, so that the choosing the initial point 
is a real problem. 

In our explicit computations, for higher periodic points (of periodicity m > 100), the 
domain of attraction (convergence) reduces so much that it appears very difficult to localize 
it (as difficult as to find the periodic point itself). 

There also exist numerical methods which converge exponentially but which have a much 
larger basin of attraction than the Newton method: minimizing methods for instance. 

9.1.1 Minimizing methods 

The initial problem, that of solving the system ( |125| ) , is equivalent to minimizing the function 
/ defined by 



The weights Wj are positive and are chosen according to dimensionality and relative precision 
criteria. In our case we simply choose u>i = 1 . 

Like in the case of solving nonlinear, non-algebraic equations, there exist only iterative 
minimization methods, which find the relative minimum closest to the starting point. The 
determination of an appropriate starting point will be a separate problem. 

The most efficient minimization methods make use as much as possible of the differen- 
tiability properties of the functions to be minimized (called the "objective function"). 

We will consider only this class of methods. In all of this class of methods, the important 
step is considered to be the choice of a direction, along which the one-dimensional minimiza- 
tion is performed, starting from a previous approximation. How to make the most efficient 
one directional minimization, is usually not treated in the literature. In contrast, we will 
treat these two problem together. 

Let d = (di, dff) be a direction, and let a = (ai, on) be some approximation of the 
minimum point, to be improved. Let us consider the following function of a single variable 
t: 



N 




(126) 



G(*)=/(a+d.t) 



(127) 



where / is given by (126) and 



x = a + d.t 



(128) 



The derivative at t — is: 



\df(*)] 



<d,(Vf) x=a > 



(129) 



dXj 
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9.1.2 The Laplace method 



The oldest method is the Laplace, or the so-called steepest descent method. The Laplace 
method consists in the following choice 

d = -(V/) x=a (130) 

such that G'(t = 0) < and that G decreases along positive values of t, at least at the 
beginning. A Laplace iteration consists in finding the minimum t\ along this direction, in 
moving to this new point ai= a + d.ti, in computing a new direction di= — (V/) x=ai by 



(130) in the new point and in repeating this procedure up to desired precision is reached. 
Clearly this method cannot be the best one because it is not invariant with respect to 



the coordinate changes, even for linear transforms: from (127) d must be a contravariant 
vector, while from ( |130| ) it is a covariant one. This lack of invariance appears more clearly in 
the case where the Xi variable components are physical quantities of different kind, so that 
this Laplace method cannot be optimal. Nevertheless, the method works exactly, if /(x) is 
of the form 

f(x.)=gi(x 1 ) + ..+g N (x N ) (131) 
where the gi are functions of a single variable. 



Another drawback of the Laplace method is related with the fact that even in the case 
where /(x) is a polynomial of second order, it remains approximative, i.e. it requires an 
infinite number of iterations to reach the exact result; excepted in the case where it is of the 



form (131). Even worse, the rate of convergence is of the order 0((A — A)/(A + A)) where 
A and A are the largest and smallest eigenvalues of the associated quadratic form (they 
are nonnegative, otherwise the minimum would be — oo ). If A>A the convergence rate is 
very slow. 



The above mentioned problems can be solved by another class of methods, known as 
conjugated direction method" ( the etymology will be explained below). 



9.1.3 Conjugated direction method 

The methods in this class have the following properties: 

A) If /(x) is a second order polynomial, then the method is exact in N steps . 

B) For computing the directions of successive one-dimensional minimizations, at most 
the first order derivatives are used. 



An important consequence of property A) is the very quick rate of convergence, in the 
general, nondegenerate case: the error decreases like 

0(exp(-fc.2 p )) (132) 

where k is a positive constant depending on initial guess and p.N is the total number of 
iterations. 



This conjugated direction method consists in : 

(i) approximating locally, at each iteration step, the objective function by a second order 
polynomial 

fii) finding a linear change of variables so that this polynomial becomes of the form of 



79 



The step (ii) is equivalent to finding a new coordinate system, generated by the conju- 
gated directions di,d 2 , ...,d N obtained by the ./V successive minimizations, such that 

(di .Q.d,-) =0 if i£j (133) 

where Q is the Hessian matrix, i.e. is the matrix of the second order derivatives of the 
objective function. The relation ( |133| ) can be reinterpreted in the following way. The Taylor 
approximation of the objective function to be minimized can be written as 



b ,h> 

d 2 f 



/(x +h) = /(x ) 

where b =V/(x ) = (V/) Xq and Q : 
defines a new scalar product denoted by square brackets 

[x,y] = (x ,Q.y) 



(h,Q.h>/2 +0(|/i| £ 



(134) 



dxidxj 



In Eq.(134) the quadratic term (h , Q.h) 

(135) 



The relation (133) expresses the orthogonality [di,dL] = of the conjugated directions 
obtained at the i — th and j — th minimization, relative to this new scalar product. Geomet- 
rically this means that in the new coordinate system, defined by the conjugated directions 
(see below how to generate these ones), the level surfaces (let us say in 3 dimensions) are 
homothetic ellipsoids, with the same center, proportional semi-axis, and axes along coordi- 
nate axes. The conjugated directions are the coordinate axes. Clearly, the minimum (the 
center) is reached by a succession of 3 steps, minimizing along coordinates. 

From the above explanation it follows that the 2 d order Taylor expansion term in Eq. 



( 134 ) of the objective function generates (locally) a natural, canonical Euclidean geometry, 
in the space of variables which can be of very different significance. 

Remark: Unfortunately, all of these general purpose algorithms are local. They require 
an initial point to start in order to find some local minimum in the neighborhood ( in fact 
a basin of attraction or " hydrographic basin" to which the initial point belongs; their 
boundaries are regular, not fractal like the basin of convergence of the Newton iterations). 



9.1.4 The Fletcher- Reeves method 

There are different ways to generate the conjugated directions. For most of applications 
the optimal one is the so called Fletcher- Reeves method ]78| . The optimality refers to 
the reduced number of operations and storage requirement of order O(N) to generate the 
conjugated directions, by a recursive process, using the gradients of the objective function. 



One main iterative step (which minimizes exactly a quadratic function) consists in N 
sub-steps. Each sub- step, indexed by k, denoted by Sk consists in the following operations. 
Denote by x the point obtained in previous step or just the initial guess of the minimum 
point, by do, the first conjugated direction. 

Sub-step So : At the first substep (fc = ) we take 

d = -V/(x ) (136) 

Compute the point X! by minimizing the function /(x + t. d ) = ipo(t) with respect to 
variable t. This is performed here by using simply the Newton method (see next Section). 
Let the optimal value be t* . Then xi = xq + t* . do 

Suppose we performed the sub-steps Sq , Sk, let dfc the direction computed in sub-step 
Sk and let x^+i be the point reached at the end of sub-step Sk- 
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Sub-step Sk+i (k < N) In the point reached in the previous sub-step, a new conjugated 
direction, dfc+i, is computed, according to: 



d fc +i = -V/(x fc+1 )+/3 fc . d fc (137) 

where (3 k is defined by 



= ||V/(x fc+1 ,H 
l|V/(x fc )|| 2 

in which the norm is the standard Euclidean one. 

Then we find the new point x^+2 by minimizing the function 

/(x fc+ i+t.d fc+ i)=V fc+ i(*) ( 139 ) 
with respect to the t variable, and using the optimal value t* . Then we compute n.k+2 by 

x/c+2 = x fe+ i + t*. d k+1 (140) 

These sub-steps are continued recursively, while the sub-steps Sq , Si,... , Sjv-i ar< 3 
executed. After performing the sub-step Sn-i the iteration step is completed and the 
approximate minimum point xjv is reached. If the objective function is a quadratic one, 
the point xjv is the exact minimum point, reached just in ONE iteration step, with iV 
sub-steps. If the objective function is not quadratic, then we go to the next iteration step, 
by using for starting point at sub-step So the last obtained point xjv, until the required 
precision is reached. 

It was prove d in J |7q that the so-computed directions are really conjugated to each other, 
i.e. they satisfy (|l33|). 
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We have used this method to find periodic points, but by applying moreover some 
improvement in the most frequent (consequently the most time-consuming) steps like the 
one-dimensional minimization. 



9.1.5 The necessity to improve the one-dimensional minimization 



By Eq.( |l32| ) - where p is the number of full iteration steps (each step involving N 
substeps)-, the speed of convergence is very high, let us say like a double exponential 
(exponential of exponential in the number of iterations). In contrast, if we use for the one- 
dimensional minimization only the values of the function (/(x& + t. d&) = tp(t)) and of the 
first order derivatives, the speed of convergence in the one-directional Newton minimization 
is only (simply) exponential. Consequently the overall convergence will be slowed down to 
simple exponential rate. 

The general minimization problem for the function /(x^ +t. dk) = V^M Q139Q can be 
reduced to solving the equation 

v4W = o (ui) 

and by Newton method one obtains ip' k (ti) + ijj'l{ti) .{ti+\ —ti) — which requires to compute 
the second order partial derivatives of the objective function. If we can compute ip" (t), then 



81 



the one-dimensional optimization can be accelerated by the one-dimensional Newton method 
which gives: 



ti 



U 



v4(*0/<(*0 



(142) 



where the index i labels the Newton iterations at substep k. In the neighborhood of an 
optimal point its convergence is very fast, i.e double exponential with respect to the number 
of Newton iterations. 

Nevertheless in our problem, concerning high order periodic points, this method is very 
difficult to use: in order to compute ip (tfc) we must compute the second order derivatives 
of the m— times iterated map T (the Tokamap for instance [H), which leads not only to 
a complicated program, but rather to a very CPU time-consuming program in the case of 
large values of m. 

9.2 The improved one-dimensional minimization 

We will see how to avoid the computation of second order derivatives, without destroying 
the double exponential convergence. This is the main improvement we bring to the Fletcher- 
Reeves method. Such a "trick" works only in the case of the minimization problems arising 
from the reduction method, discussed above: we will suppose that the functional to be 
minimized is of the form of Eq. (126). Let ip(t) be one of the functions ^fc(i), and d one of 
the directions d^. 



*W=/(x + t. d) 
where the vector d is one of the conjugated directions. Then 



OF* 



dxi 



.Fj(x + t. d).di 



(143) 



(144) 



x+t.d 



and 



where 



d 2 i/j/dt 2 = 2 



g(t) 



ait) 



N 

E ■ 



N 

E 

i,j,fc=i 



d 2 F 3 
dxi.dxk 



dxi dxk 



■ di di 



-id 



.F Q (x + t.d).di d k 



-t.d 



(145) 



(146) 



For the case of searching periodic points, the computation of first order derivatives reduces 
to the computation of the products of Jacobi matrices. On the other hand, the second 



term (146) Q x .g 3 Xk -Fj is small near the solution, because of its last factor, and will be ne- 
glected. It can be proved that the so-obtained iterative process still has a double exponential 
convergence. 



9.2.1 Summary of the improved Fletcher Reeves method. 

Shortly, the essence of our method is the following: The optimal value of t* where the 
minimum of 



/(x fe+ i + 1. dfc+i) = ip k+1 (t) = ip(t) 



(147) 
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is reached, is approximated by 



V>'(0) 



(148) 



A(0) 



where i/j (t) is given by ( 144 ) and where we have introduced 



N 




m = E 



(149) 



where the dk+i,i are defined from the I th component (I = 1,2, ...,N) of the vector at step 
k + 1, i.e. from: 



In the rest, the steps follows the standard Fletcher-Reeves strategy. Our modification: 

(a) preserves the same double exponential convergence speed and 

(b) saves computer CPU time, for the following reasons: 

(bl) we do not compute the second derivatives of Fj. These computations in the 
case of periodic point search require CPU times which increase quadratically with m, the 
periodicity order 

(b2) we perform SINGLE iterations, instead of many (at least four) Newton iterations 
at each substep. The CPU time is reduced at least by a factor of four. 

(c) is not unstable, do not produce overflow and has larger domain of convergence, 
compared to Newton method. 

We observed a single drawback: the tendency to remain "suspended" in some local 
minima, for larger order periodic points (hundreds). 

This problem was corrected by refining the output (i.e. periodic point coordinates) of the 
conjugated gradient program by another new program. This program now uses the classical 
Newton method. The instability of the Newton method, mentioned above, does not show 
up, because the starting point (initial guess for the Newton method) is already very close 
to the solution and it is in the domain of convergence (or attraction basin). 

9.2.2 General strategy for periodic point search (programming details) 

In order to localize a periodic points of periodicity m and having a given type (n, m) we 
must realize the following steps: 

1) to find an initial point for the minimizing method. 

This is achieved in our program by a preliminary search for the absolute minimum of the 
objective function on a finite grid inside the input tetragon. The coordinates of the grid's 
point giving the lowest value serve as input for conjugated gradient search. 

The vertices of this tetragon, for the lowest order periodic points, are obtained from 
preliminary graphics of tokamap phase portrait. For highest order, the input tetragon was 
obtained recurrently form previous lowest order periodic points. 

2) to minimize the objective function, using the improved Fletcher-Reeves method, with 
the starting point determined by 1). 

For this we used a C-\ — |- program which contains the improved Fletcher-Reeves method. 

If the output error of the conjugate gradient program is too large (as compared with a 
maximal imposed error), the results can be refined by the Newton iteration program, the 
input of which is the output of previous program. 

The output of the program contain: the coordinates of the searched periodic point. 



(dk+i,i, — i dk+i.i, dk+i.N) 



(150) 
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The program computes the safety factor related to the magnetic axis (i.e. m and n) 
and the error (by comparing the initial and final positions of the chain of periodic points). 
The cross-check of numerical floating-point accuracy is the deviation of the value of the 
determinant of the Jacobi matrix from its known unity value (in the output it is named as 
"determinant" error). The complete program files can be obtained by sending a mail to one 
of the authors []. 

The first program computes the safety factor related to the magnetic axis, i.e. m and n. 
This program, with simple modifications, can be used in different cases, in order: 

a) to find periodic points of 2 dimensional maps, 

b) to solve general nonlinear system of equations, 

c) in case of field- line tracking, to find the periodic field lines. In this case the differential 
equations of the field lines must be completed with the so called Jacobi variational equations 
(i.e. the linearized equations), or more generally 

d) to find the periodic trajectories of the dynamical systems with finite degree of freedom, 
described by discrete iterative maps or by differential equations. 
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